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Abstract 

A nontraditional numerical method for solving 
conservation laws is being developed. The new method 
is designed from a physicist’s perspective, i.e., its de- 
velopment is based more on physics than numerics. 
Even though it uses only the simplest approximation 
techniques, a 2D time-marching Euler solver developed 
recently using the new method is capable of generating 
nearly perfect solutions for a 2D shock reflection prob- 
lem used by Helen Yee and others. Moreover, a recent 
application of this solver to computational aeroacous- 
tics (CAA) problems reveals that: (i) accuracy of its 
results is comparable to that of a 6th order compact 
difference scheme even though nominally the current 
solver is only of 2nd-order accuracy; (ii) generally, the 
non-reflecting boundary condition can be implemented 
in a simple way without involving characteristic vari- 
ables; and (iii) most importantly, the current solver is 
capable of handling both continuous and discontinuous 
flows very well and thus provides a unique numerical 
tool for solving those flow problems where the interac- 
tions between sound waves and shocks are important, 
such as the noise field around a supersonic over- or 
under-expansion jet. 
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Figure 1 . — Pressure distribution for a shock reflection 
problem (a flow of Mach number 2.9 enters from the left; 
the shock is reflected by a wall on the back). 


1. Introduction 

The method of space-time conservation element 
and solution element is a new numerical discretiza- 
tion method for solving conservation laws [1-13]. The 
new method differs substantially in both concept and 
methodology from the well-established methods-i.e., 
finite difference, finite volume, finite element, and 
spectral methods. It is designed to overcome sev- 
eral key limitations of the above traditional methods. 
In addition, its development is driven by an empha- 
sis on simplicity, generality, and accuracy. In this 
paper, we shall describe a 2D time-marching Euler 
solver developed recently using the new method [2]. 
Even though it does not use (i) any approximation 
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techniques more complicated than Taylor’s expansion 
or weighted averaging, (ii) any mesh-refinement tech- 
niques, (iii) any monotonicity constraints, (iv) any 
characteristics-based techniques, or (v) any ad hoc 
techniques that are used only in the neighborhood 
of a discontinuity, this solver is capable of generat- 
ing highly accurate solutions for a 2D shock reflection 
problem used by Helen Yee and others [14]. As shown 
in Fig. 1, both the incident and the reflected shocks can 
be resolved by a single data point without the presence 
of numerical oscillations near the discontinuity . 

Recently, Loh and others [13] have successfully 
simulated some benchmark aeroacoustics problems us- 
ing the new Euler solver. The non-reflecting boundary 
condition is nearly perfect with almost no modification 
of the original code. In several free shear layer calcu- 
lations, the smallest eddies are consistently resolved 
in 3-4 grid cells. Some interesting vortex/shock (or 
Mach waves) interactions, which have never been re- 
ported before, were also observed. A further discussion 
of these results is given near the end of Sec. 5. 

The Introduction section of [1] begins with a 
lengthy discussion that includes (i) an examination of 
the traditional methods and their limitations, and (ii) 
a study of the general requirements for an accurate 
solver of conservation laws. The current method is 
built on a set of design principles that emerge from 
the above discussion. In this section, an overall view 
of the current method will be given from a historical 
perspective. Particularily, it will be explained how the 
current method was shaped by fundamental physics 
considerations. 

The first physical problem considered in the cur- 
rent development is an initial- value problem involving 


the PDE 


du du 

» +<, 9 ? = 0 


( 11 ) 


where the convection speed a is a constant. The exact 
solution to any such problem has three fundamental 
properties: (i) it does not dissipate with time, (ii) its 
value at a spatial point at a later time has a finite do- 
main of dependence (a point) at an earlier time, and 
(iii) it is completely determined by the initial data at a 
given time. Ideally, a numerical solution for Eq. (1.1) 
should also possess the same three properties. Because 
(i) a solution of a dissipative numerical scheme will 
dissipate with time; (ii) the value of a solution of an 
implicit scheme at any point (a?,f) is dependent on all 
initial data, and all the boundary data up to the time 
t ; and (iii) the unique determination of a solution of 
a multi-level scheme requires the specification of the 
initial data at two or more time levels, an ideal solver 
must be a two-level, explicit, and neutrally stable (i.e., 


non-dissipative) scheme . In 1991, such a solver was re- 
ported in [5]. Because this new solver models Eq. (1.1) 
which is characterized by the parameter a, hereafter it 
will be referred to as the a scheme. The a scheme is 
neutrally stable if the Courant number < 1. It is the 
only two-level explicit solver of Eq. (1.1) known to the 
authors to be neutrally stable. 

Development of the a scheme is aided by recasting 
Eq. (1.1) into an integral form that represents a law of 
flux conservation [1,5,9]. The a scheme is derived by 
using both Eq. (1.1) and its integral form. Note that 
the conservation law appears in a form in which space 
and time are unified and treated on the the same foot- 
ing. This unity of space and time, and the requirement 
that Bux-conservation be enforced locally (i.e., down to 
a computational cell) and globally (i.e., over the entire 
computational domain) are two tenets of the current 
numerical method. They are also the key characteris- 
tics that distinguish the current method from most of 
the traditional methods. 

Note that the enforcement of both local and 
global Bux conservation may also prevent inconsis- 
tency among initial/boundary conditions. 

The a scheme has many nontraditional features. 
At each mesh point ( j,n ), it has two independent 
marching variables ti? and (ti r )j with the latter be- 
ing the analogue of du/dx at the mesh point. It also 
has the simplest stencil, i.e., a triangle with one ver- 
tex at the given time level and the other two vertices 
at the previous time level. In contrast to a typical 
finite- volume scheme, extrapolation and interpolation 
are not used in the a scheme for the evaluation of the 
flux at an interface separating two neighboring conser- 
vation elements. Moreover, the a scheme is a two-way 
marching scheme, i.e., the forward marching scheme 
can be inverted to become the backward marching 
scheme. In other words, the marching variables at 
the (n — l)th time level can be determined in terms 
of those at the nth time level. These and other non- 
traditional features of the a scheme are discussed in 
depth in [1,9]. 

Because there are two independent marching vari- 
ables at each mesh point (j, n), two amplification fac- 
tors appear in the von Neumann stability analysis of 
the a scheme [1,5,9]. It happens that these two factors 
are identical to those of the Leapfrog scheme [15] if the 
latter factors arise from a “proper” von Neumann anal- 
ysis. Note that a solution to the main Leapfrog scheme 
(excluding the starting scheme which relates the mesh 
variables at the first two time levels) is formed by two 
decoupled solutions. Traditionally the von Neumann 
analysis for the Leapfrog scheme is performed without 
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taking into account this decoupled nature. In [1,9], 
it is performed separately for each decoupled solution. 
The amplification factors thus obtained are identical 
to those of the a scheme. This coincidence was un- 
expected because the Leapfrog and the a schemes are 
structually different. 

The fact that the amplification factors of the 
a scheme are related to those of a celebrated clas- 
sical scheme is only one among a string of similar 
unexpected coincidences encountered during the de- 
velopment of the current method . As will be ex- 
plained in this paper , the amplification factors of the 
Lax-Wendroff t the Lax , the Crank-Nicolson , and the 
DuFort-Frankel schemes [15] also are identical to those 
of some of the extensions of the a scheme. 

Because the solutions of the a scheme were re- 
quired to share with those of Eq. (1.1) the same funda- 
mental properties referred to earlier, during the early 
day of its development Chang realized that such a 
scheme must share with Eq. (1.1) certain basic in- 
variant properties in space- time. He studied the in- 
variant properties of Eq. (1.1) with respect to spa- 
tial reflection, time reversal, and space-time inversion. 
He then defined rigorously what it means for a two- 
level constant-coefficient finite-difference analogue of 
Eq. (1.1) to have similar invariant properties. In this 
study, it is assumed that there is only one marching 
variable ti* at a mesh point for a finite-difference ana- 
logue of Eq. (1.1). With the above assumptions, it 
was shown that the von Neumann amplification factor 
of a numerical analogue satisfies a special relation for 
each invariant property that this analogue possesses. 
Particularly, an analogue is unconditionally neutrally 
stable if it is invariant under space-time inversion. As 
an example, the implicit Wendroff scheme [16] is in- 
variant under space-time inversion. Its amplification 
factor is 


_ cos(fl/2) - tVsm(g/2) _ / — 
~ cos(0/2) + iv sin(0/2) V V ’ 


( 1 . 2 ) 


Here 6 is the phase angle variation per mesh interval 
and u is the Courant number. Note that |G| = 1 for 
any v. 

The results of the above study were presented in 
a conference that took place in 1992 [6]. Because a 
two-level explicit finite-difference scheme is not invari- 
ant under space- time inversion, one can infer from the 
arguments made in [6] that such a scheme cannot be 
an ideal scheme, i.e., neutrally stable. 

The two-level finite-difference scheme has only 
one amplification factor. For a numerical analogue of 
Eq. (1.1) with two amplification factors G+ and G_, 


generally one cannot conclude that the scheme is un- 
conditionally neutrally stable if it is invariant under 
space-time inversion. For the Leapfrog scheme or the 
a scheme, we have G+G- = —1. In general it does not 
imply neutral stability, i.e., |G+| = |G_| = 1. How- 
ever, it does imply that such a scheme must be neu- 
trally stable if it is stable , i.e., if |G+| < 1 and |G_ | < 
1. For both the a and the Leapfrog schemes, this im- 
plies that they are neutrally stable if the Courant num- 
ber < 1. 

In the 11th AIAA CFD conference (July, 1993), 
Thomas and Roe* presented a paper [17] in which 
the concept of the invariance under space-time inver- 
sion is also used to construct non-dissipative numerical 
schemes. Note that the concept of “rotational sym- 
metry” in space-time (i.e., “invariance under space- 
time inversion”) was discussed briefly in [17]. Unfor- 
tunately, this discussion is faulty in two aspects. It 
considers only one amplification factor G in spite of the 
fact that the scheme under consideration, i.e., the “Up- 
wind Leapfrog scheme” has two amplification factors. 
Furthermore, it leads to the conclusion that G = 1 for 
a scheme that is invariant under space-time inversion. 
Even for a scheme that has only one amplification fac- 
tor, this conclusion is obviously false. As an example, 
consider the Wendroff scheme. It is invariant under 
space-time inversion and yet G ^ 1, except for the 
case v sin(0/2) = 0 (see Eq. (1.2)). 

In a 1994 ICASE report [18], the concept of in- 
variance under space-time inversion is again used by 
Roe to construct non-dissipative linear Bicharacteris- 
tic schemes. 

The a scheme is only a special case of a more 
general two-level explicit scheme described in [5]. It is 
a solver for 


du du d 7 u 

dt +a fo- ti jh* 


= 0 


(1.3) 


where the viscosity coefficient p(> 0) is a constant. 
Because this solver models Eq. (1.3) which is charac- 
terized by the parameters a and /i, hereafter it is re- 
ferred to as the a-p scheme. The a-p scheme reduces 
to the a scheme if p = 0. Because the a scheme is neu- 
trally stable, the a-p scheme has the property that the 
numerical dissipation of its solution approaches zero 
as the physical dissipation approaches zero. 

*A participant of the 23rd Conference on Modeling 
and Simulation, April 30-May 1, 1992, Pittsburgh, PA. 
He was in the audience when the paper [6] was pre- 
sented by Chang, and participated in the discussion 
afterward. 
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The above property is important because of the 
following observation: With a few exceptions, numer- 
ical dissipation generally appears in a numerical so- 
lution of a time-marching problem. In other words, 
the numerical solution dissipates faster than the corre- 
sponding physical solution. For a nearly inviscid prob- 
lem, e.g., flow at a large Reynolds number, this could 
be a serious difficulty because numerical dissipation 
may overwhelm physical dissipation and cause a com- 
plete distortion of solutions. To avoid such a difficulty, 
at the minimum, a model solver for Eq. (1.3) should be 
required to have the special property mentioned above. 

In [5], it is shown that the a-/i scheme, which is 
explicit, has the unusual property that its stability is 
limited only by the CFL condition, i.e., it is indepen- 
dent of /i. Furthermore, the. amplification factors of 
the a-p scheme reduce to those of the DuFort-Frankel 
scheme if a = 0 [1,9]. Note that a solution of the latter 
scheme is also formed by two decoupled solutions. As a 
result, the von Neumann analysis should be performed 
separately for each decoupled solution. 

In order to “explore the concept of a dynamic 
space-time mesh and the need for a unified treatment 
of physical variables and mesh parameters” [5, p.24], 
a moving mesh with a uniform speed b was introduced 
in [5] (see Fig. 2(a)). With the aid of a Galilean trans- 
formation, it is shown that the a-p scheme or a typi- 
cal finite-difference solver of Eq. (1.1) or (1.3) can be 
converted to its moving-mesh form by simply replac- 
ing the parameter a with a — fc. As an example, the 
moving-mesh form of the Leapfrog scheme is 


..n+i _ t #”" 1 *# n — .,n 

^r x - + (°-«- i± k ^= 0 


Note that: (i) Eq. (1.4) reduces to the regular form 
of the Leapfrog mesh if 6 = 0, i.e., if the mesh be- 
comes stationary, and (ii) a moving-mesh form and 
its stationary-mesh counterpart actually represent two 
different schemes (read Sec. 3, “The Dynamic Space- 
Time Mesh” in [5]). 

Let b = a. Then (i) a mesh line with j being a 
constant along this line points in the characteristic di- 
rection of Eq. (1.1), and (ii) Eq. (1.4) is reduced to 
ti” +1 = tij” 1 , i.e., the numerical value of u is con- 
stant along such a mesh line. Because u is constant 
along a characteristic line if it is an exact solution of 
Eq. (1.1), aside from round-off errors, a numerical solu- 
tion matches perfectly with its analytical counterpart. 

As a result of the above observation, and other 
considerations given in [5], one arrives at the impor- 
tant conclusion [5, p.3]: “• • • (i) stability and accuracy 
can be improved , and (ii) dissipation and dispersion 


can be reduced , if the space-time mesh is allowed to 
evolve with the physical variables such that the local 
convective motion of physical variables relative to the 
moving mesh is kept to a minimum .” 

Because of the complexity involved, the dynamic 
space-time mesh has not yet been used in the current 
development beyond the applications reported in [5]. 

Note that the “Upwind Leapfrog scheme” referred 
to in [17,18] can be considered as a special case of 
Eq. (1.4). Let b = A x/a*. Then, as depicted in 
Figs. 2(a) and 2(b), the mesh points (j,n + 1) and 
(j + l,n) line up at one spatial location, while the 
mesh points (j — 1, n) and (.;, n — 1) line up at another 
spatial position. Also Eq. (1.4) becomes 

w * 1 --?*■> +(»?-i -«r‘) . 

2a* 2ax 

(1.5) 

The stencil of Eq. (1.5) is shown in Fig. 2(b). The 
stencil of the “Upwind Leapfrog scheme” is depicted 
in Fig. 4 in [17] and reproduced in Fig. 2(c) here. 

A comparison between Figs. 2(b) and 2(c) reveals 
that (i) the mesh points (j, n 4- 1), (j -f 1, n), (,; — 1, n), 
and (j,n — 1) in Fig. 2(b) correspond, respectively, 
to the mesh points (j, n -f 1), (j,n), ( j — l,n), and 
(j — l,n — 1) in Fig. 2(c); (ii) the mesh interval 2ax 
in Fig. 2(b) corresponds to the mesh interval ax in 
Fig. 2(c); and (iii) the time-step size a t in Fig. 2(b) 
corresponds to the time-step size At in Fig. 2(c). Note 
that the moving-mesh indices are used here, while the 
stationary-mesh indices are used in [17]. By replac- 
ing (i) the moving-mesh indices with the correspond- 
ing stationary-mesh indices, and (ii) 2ax with ax, 
Eq. (1.5) becomes 


K +1 ~ “?) + K-i - ^i 1 ) , - “i-i = 0 

2a t AX 

( 1 . 6 ) 

i.e., the defining equation for the “Upwind Leapfrog 
scheme” (see Eq. (4) in [17]). 

The a scheme is neutrally stable and reversible 
in time [1,5,9]. It is well known that a neutrally sta- 
ble numerical analogue of Eq. (1.1) generally becomes 
unstable or highly dispersive when it is extended to 
model the Euler equations. It is also obvious that a 
scheme that is reversible in time cannot model a physi- 
cal problem that is irreversible in time, e.g., an inviscid 
flow problem involving shocks. Hence, the a scheme 
is extended to become the a-e scheme [1,9]. Stability 
of this scheme is limited by the CFL condition and 
0 < c < 1 where e is a special parameter that con- 
trols numerical dissipation. Moreover, if e = 0, the 
a-e scheme reduces to the a scheme which has no nu- 
merical dissipation. On the other hand, if e = 1, the 
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two amplification factors of the a-e scheme become the 
same function of the Courant number and the phase 
angle. Unexpectedly, this function also is the amplifi- 
cation factor of the highly diffusive Lax scheme. 

The a-e scheme was extended [1,9] to solve the 
ID time-dependent Euler equations of a perfect gas. 
The Euler extension has stability conditions similar to 
those of the a-e scheme itself. It has the unusual prop- 
erty that numerical dissipation at all mesh points can 
be controlled by a set of local parameters. Moreover, it 
is capable of generating accurate shock tube solutions 
with the CFL number ranging from 1.0, to 0.022. 

The a-p scheme was also extended [1,9] to solve 
the ID time-dependent Navier-Stokes equations of a 
perfect gas. Stability of this explicit solver also is lim- 
ited only by the CFL condition. In spite of the fact 
that it does not use (i) any techniques related to the 
high-resolution upwind methods, and (ii) any ad hoc 
parameter, the Navier-Stokes extension is capable of 
generating highly accurate shock tube solutions. Par- 
ticularly, for high-Reynolds-number flows, shock dis- 
continuities can be resolved within one mesh interval. 

This concludes the review of the development of 
the current method up to 1993. In the rest of this 
paper, we shall describe several 2D extensions of the 
ID schemes discussed earlier. The first scheme to be 
discussed is the 2D version of the a scheme. Again, it 
is a two-level, explicit, and neutrally stable scheme. 


2. The 2D a Scheme 

In this section, we consider a dimensionless form 
of the 2D convection equation, i.e., 


du du du 

« + °' 5J + 0 'S 


= 0 


( 2 . 1 ) 


where a r , and a y are constants. Let x\ = z, z 2 = y, 
and a? 3 = t be the coordinates of a three-dimensional 
Euclidean space E 3 . By using Gauss’ divergence the- 
orem in the space-time £3, it can be shown that 
Eq. (2.1) is the differential form of the integral con- 
servation law 

(f hds = 0 (2.2) 

Js(v ) 

Here (i) S(V) is the boundary of an arbitrary space- 
time region V in £3, (ii) 


h = (a r ti, a y ti, u) (2.3) 


is a current density vector in £3, and (iii) ds = dan 
with da and n, respectively, being the area and the 
outward unit normal of a surface element on S(V). 


Note that (i) h • ds is the space-time flux of h leav- 
ing the region V through the surface element ds } and 
(ii) all mathematical operations can be carried out 
as though £3 were an ordinary three-dimensional Eu- 
clidean space. As will be shown shortly, £3 will be di- 
vided into nonoverlapping space-time regions referred 
to as conservation elements (CEs). 

In the following, we shall introduce a nontradi- 
tional space-time mesh. Its use will result in the sim- 
plest stencil possible for the current scheme, i.e., a 
tetrahedron in 3D space-time with one vertex at the 
upper time level and other three at the lower time 
level. Also this mesh is compatible with the require- 
ment that the 2D a scheme be invariant under space- 
time inversion. Most importantly, three conservation 
elements per mesh point are embedded in this mesh 
such that each mesh point is associated with three flux- 
conservation conditions. This mesh structure matches 
perfectly with the fact that there are three unknowns 
(the numerical analogues of u, du/dx and du/dy) at 
each mesh point. 

Let n denote the time level and 

t n = riAt, n = 0, ±1/2, ±1, ±3/2, . . . (2.4) 

Let j and k be spatial mesh indices with j, k = 
0, ±1/3, ±2/3, ±1, . . . (see Figs. 3-6). Let de- 
note the set of mesh points (i, fc,n) with j, k = 
0,±1,±2, ..., and n = ±1/2, ±3/2, ±5/2,.... These 
mesh points are marked by solid circles. Let f} 2 
denote the set of mesh points (i, 4,n) with j, k = 
1/3, 1/3±1, l/3±2, . . ., and n = 0, ±1, ±2, — These 
mesh points are marked by open circles. The union of 

and fi 2 will he denoted by Q. 

Each mesh point (j, Jb,n ± 1/2) in (by defi- 
nition, this implies that n = 0,±1,±2,...) is associ- 
ated with three CEs, denoted by CE^(j,i,n + 1/2), 
£ = 1,2,3 (see Fig. 7(a)). It is also associated with a 
solution element (SE), denoted by SE ^^(i,i,n± 1/2) 
(see Fig. 7(b)). Similarly, each mesh point (j,Jfc,n+ 1) 
in Qo is associated with three conservation elements 
CE i(j,k, n ± 1), £ = 1,2, 3 (see Fig. 8(a)), and a so- 
lution element SE^(j, t,n + 1) (see Fig. 8(b)). Each 
CE is a quadrilateral cylinder in space-time while each 
SE is the union of three vertical planes, a horizontal 
plane, and their immediate neighborhood. The ge- 
ometry of the hexagon ABCDEF , which appears in 
both Figs. 7(a) and 8(a), is determined by three pos- 
itive parameters it/, b and h (see Fig. 9(a)). Without 
any loss of generality, we assume that the line segment 
joining points D and A in Fig. 9(a) is parallel to the 
z-axis. Note that the form of Eq. (2.1) will not change 
under an orthogonal transformation on the x-y plane. 
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Thus one always can introduce a set of Cartesian co- 
ordinates (x,y) such that the above line segment is 
parallel to the x-axis. However, the values of a x and 
a y may change as a result of such a transformation. 

According to Fig. 3, E 3 can be filled with the CEs 
defined above. Moreover, it is seen from Figs. 7(a), 
7(b), 8(a), and 8(b) that the boundary of a CE is 
formed by the subsets of two neighboring SEs. 

Let the space-time mesh be uniform, i.e., the par 
rameters a t, w , 6, and h are constants. Let x Jt k and 
yj t k be the x - and y- coordinates of any mesh points 
(i, Jb, n) £ Q. Let x 0> o = 0 and 3/0,0 = 0. Then infor- 
mation provided by Figs. 9(a) and 9(b) implies that 

x i,t = U + k ) w + (*“ j) b , yj.k = (*~ j) k (2.5) 


(xj, 3/*, t n ). As a result, the expression on the right side 
of Eq. (2.8) can be considered as the first-order Tay- 
lor’s expansion of u(x,y,t) at ( Xj,yt,i n ). Also note 
that Eq. (2.9) is the numerical analogue of Eq. (2.3). 

We shall require that u = u*(x,y,t\j, i,n) satis- 
fies Eq. (2.1) within SE(j,ib,n). As a result, 

(«<)", k = - + a v( u v)>,i] (2-10) 

Substituting Eq. (2.10) into Eq. (2.8), one has 

+ Mh tt* - - ««(< - *")] (2.11) 

+ («»)",* [(y - yj.k) - - HI 


Let n i, n 2, f?3, n 4, ns, and «6 be the vectors depicted 
in Fig. 9(a). They lie on the x-y plane and are the 
outward unit normals to AB, BC , CD , DE , EF , and 
FA, respectively. It can be shown that 


f? - (ft> —6 + W3, 0) 
7,1 _ y/h*+(b-w/ 3) 2 ’ 


S4 = — r»j (2.6a) 


Thus there are three independent marching variables, 
i.e., u” t , (u*)" t , and (uj,)" t associated with a mesh 
point ( j,k,n ) € ft. For any (j,k,n + 1/2) € fti, 
these variables will be determined in terms of those 
associated with the mesh points (j + 1/3, Jt + 1/3, n), 
(J - 2/3, it + 1/3, n), and (j + 1/3, k — 2/3, n) (see 
Fig. 10(a)) by using the flux conservation relations: 


n 2 = (0,1,0), rj 5 = -n 2 


and 


( 2 . 66 ) 


<f h'-ds = 0, ^ = 1,2,3 (2.12) 

Js(CE[ l) (j,k, n+1/2)) 


«3 = 


(— /», 6 4- tn/3,0) 

Vh 2 + (6 + u)/3) 2 ’ 


n 6 = — n 3 (2.6c) 


For any (j, Jb, n) € ft, let 

( k, n), if (j, k, n) € ft x 

SE(j, k, n) = f < 

,SE ( 2 )(i,fc,n), if (j,k,n) € ft 2 . 

(2.7) 

For any (r,y,t) 6 SE(j,fc,n), u(r,y,t) and h(x,y,t), 
respectively, are approximated by 


v'(x,y,i;j,k,n) = + («x)” lt (a: - *y.Jt) , 2 

+ (“y)",t(y - yj.t) + («<)”,*(< - O 

and 

ft 2/> £* J) *0 = (®i I/i ^5 j\ ft) *0, (2 9) 

a y u*{x, y, <; j, *, n), u*(x, y, t; j, k, n)] 

where (t/ r )" >Jk , ( w y)j t jki an< ^ (ut)j,fc ^ constants 
within SE(j, Jfc,n). The last four coefficients, respec- 
tively, can be considered as the numerical analogues 
of the values of u, du/dx , du/dy, and dufdt at 


Similarly, the marching variables at any (j,k, n + 
1) £ SI 2 are determined in terms of those associ- 
ated with the mesh points (j — 1/3, ib — 1/3, n + 1/2), 
(i+2/3,4-1/3, n+l/2), and (j- 1/3, *+2/3, n+l/2) 
(see Fig. 10(b)) by using the flux conservation rela- 
tions: 

/ h*-ds = 0, £ = 1,2,3 (2.13) 

JS(CE< a) (j,k,n+l)) 

Obviously, Eqs. (2.12) and (2.13) are the numerical 
analogues of Eq. (2.2). 

As a result of Eqs. (2.12) and (2.13), the total 
flux leaving the boundary of any CE is zero. Because 
the flux at any interface separating two neighboring 
CEs is calculated using the information from a single 
SE, the flux entering one of these CEs is equal to that 
leaving another. It follows that the local conservation 
conditions Eqs. (2.12) and (2.13) will lead to a global 
conservation condition, i.e., the total flux leaving the 
boundary of any space-time region that is the union of 
any combination of CEs will also vanish. 

In the following, several preliminaries will be given 
prior to the evaluation of Eqs. (2.12) and (2.13). To 
proceed, note that a mesh line with j and n being con- 
stant or a mesh line with k and n being constant is not 
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aligned with the x-axis or the y-axis. We shall intro- 
duce a new spatial coordinate system (£, 77 ) with its 
axes aligned with the above mesh lines (see Fig. 9(c)). 

Let e x and e y be the unit vectors in the x- and the 
y- directions, respectively. Let and t ' n be the unit 

vectors in the directions of DF and DB (i.e., the j- 
and the k- directions-see Figs. 9(a)-(c)), respectively. 
It can be shown that 



= [(^ “ 6 )e* - he y ] /a( 

(2.14) 

and 

where 

e n = [(u> + b)e x + he y ] /at) 

(2.15) 


With the aid of Eqs. (2.5), (2.20), and (2.22), it 
can be shown that the coordinates (£, 77 ) of any mesh 
point (j, fc, n) G are given by 


C = j aC, and ry = k Ary (2.23) 

i.e., aC and at; are the mesh intervals in the £- and 
the 77 - directions, respectively. 

Next we shall introduce several coefficients that 
are tied to the coordinate system (£, 77 ). Let 

(:)* 

Also, for any (j, ik,n) G ft, let 



AC d = \dP\ = V(u; - 6) 2 + h 2 


and 


= f \DB\ = y/(w + 6) 2 + h 2 


(2.16) 

(2.17) 


Let the origin of (a;, y) also be that of 77). Then the 
spatial coordinates (ac, y) and (C, tj) of any point in £3 
are related by the condition 


C e ( + 7 e v = x e z + y e v 


(2.18) 


Substituting Eqs. (2.14) and (2.15) into Eq. (2.18), one 
has 

(2.19) 

( 2 . 20 ) 


and 


Here 


and 


(;)-($) 


T = f 


,)= T -' 

(l) 

/ w — 6 

w + b\ 

A< 

A 17 

h 

h 

AC 

AT) ) 


( 2 . 21 ) 


1 dcf 


aC (w + 6)aC \ 
2u; 2wh 


Ary (tv — 6) Ary 
2wh 


' 2 u> 


( 2 . 22 ) 



where 7* is the transpose of T . For those who are 
familiar with tensor analysis, the following comments 
will clarify the meaning of the above definitions: 

(a) (a^ , ay ) are the contravariant components with re- 
spect to the coordinates (£, ry) for the spatial vec- 
tor whose x- and y - components are a x and a y , 
respectively. 

(b) ((ii( are the covariant components 
with respect to the coordinates (£ ,ry) for the spa- 
tial vector whose x- and y- components are (tir)£t 
and (%)£*, respectively. 

(c) Because the contraction of the contravariant com- 
ponents of a vector and the covariant components 
of another is a scalar, Eq. (2.10) can be rewritten 
as 


( u *)j,k = - ( u n)?,i] (2.26) 

(d) Under the linear coordinate transformation de- 
fined by Eqs. (2.19) and (2.20), (C — i aC, ry - ib a ry) 
are the contravariant components with respect to 
the coordinates (£, ry) for the spatial vector whose 
x- and y- components are x — Xj t k and y - yy f *, 
respectively. Using the same reason given in (c), 
Eq. (2.11) implies that 

u m (x,y,t;j,k,n) = U*(C, V, t;j,k,n) (2.27) 

where 

«*(<>«M;i.M) = f «",* 

+ («c)",t KC - M) - a ( (t - <")] (2.28) 

+ ( «>>)",* [(*7 - *ajj) - a„(t - < n )] 


Note that the existence of T~ l , the inverse of T, is 
assured if wh ^ 0. 
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Note that Eqs. (2.26) and (2.27) can also be ver- 
ified directly using Eqs. (2.20), (2.22), (2.24), and 
(2.25). 

Next, let (i) 


j. def 6 


(«) 


K + )",* = X <«C>jV = ( 2 “) 


and (iii) 


def At . def A/ , 

v < = T a <’ 1/17 = T a * 


(2.31) 


The coefficients defined in Eqs. (2.29) and (2.30) can 
be considered as the normsdized counterparts of those 
defined in Eqs. (2.24) and (2.25). Also note that 
and at), respectively, are the lengths of the sides DF , 
and BD of A BDF depicted in Figs. 9(a)-(c). More- 
over, by substituting Eq. (2.29) into Eq. (2.31), one 
has 



AT] 



(2.32) 


In other words, (2/3)i/{ and (2/3)^ are the Courant 
numbers in the £- and r/- directions, respectively. 
Furthermore, let be defined by 



(1)± def 1 

<T n — 1 — l<C — 

(2.33) 

4?* d = ± ( 1 - + ^c) 

(2.34) 

^3 }± = f ±(1 ~v<~ "nX 1 + "n) 

(2.35) 


(1)± def , , 
0-21 = 1 + 

(2.36) 


a 22± — T(1 + t 'c)(2 — 

(2.37) 


<*23 * = f ±(1 + V()(l + Vr,) 

(2.38) 


(1)± def - . 
<*31 = 1 + 

(2.39) 


= f ±(1 + H,)(l + V() 

(2.40) 


= ^ 1 +^x 2 -^) 

(2.41) 


<* ( u )±d =l +*'<+«, 

(2.42) 

— T(1 + + « / f ) )(l - 

(2.43) 

^3 }± = f =F(1 + "C + ^)(1 - V„) 

(2.44) 


<£ )±d = f l-c 

(2.45) 

= f ±(1 - i/ c )(2 + v ( ) 

(2.46) 

^ d ==F(l-^)(l-^) 

(2.47) 

^ )±d = rf l-^ 

(2.48) 

^32 ):t ^ =F(1 - H|)(l “ V() 

(2.49) 

<4l )±d =±( i-^)(2+^) 

(2.50) 


Note that: 

(a) Each of Eqs. (2.33)-(2.50) represents two equa- 
tions. One corresponds to the upper signs while 
the other, to the lower signs. 

(b) The definitions given in Eqs. (2.33)-(2.41) will be 

used in the first marching step of the 2D a scheme; 
while those given in Eqs. (2.42)-(2.50) will be used 
in the second marching step. It is seen that the 
expressions on the right sides of the former can 
be converted to those of the latter, respectively, 
by reversing the and a — ” signs. Moreover, 
for every pair of m and £, and are 

converted to and respectively, if i/(, 

and Vr) are replaced by — , and — i/,,, respectively. 

To simplify the following development, let 


(;, fc; 1, 1) == j + 1/3, k + 1/3 

(2.51a) 

(j, 1. 2) = f j — 2/3, k + 1/3 

(2.516) 

(;*,*; 1, 3) "; + 1/3, *-2/3 

(2.51c) 

(i, 2, l) = f j — 1/3, * — 1/3 

(2.52a) 

0", 1:; 2, 2) = f j + 2/3, k — 1/3 

(2.526) 

0, *; 2, 3) = f j — 1/3, * + 2/3 

(2.52c) 


Note that (i) (.;, Jb; 1, £), l = 1,2, 3, are the spatial mesh 
indices of points A , C, and E depicted in Fig. 7(a), 
respectively, (ii) (j,k] 2,£), l = 1,2,3, are the spa- 
tial mesh indices of points D, F , and B depicted in 
Fig. 8(a), respectively, and (iii) the mesh indices on 
the right sides of Eqs. (2.51a, b,c) cam be converted to 
those in Eqs. (2.52a, b,c) by reversing the a +” and u — v 
signs. 

Equations (2.12) and (2.13) are evaluated in [2]. 
With the aid of the above definitions, the results are 
summarized as follows: 
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(a) Eq. (2.12) with 1 = 1: 

[«8 H « 

= [°u ) ~« 

(b) Eq. (2.12) with l = 2 


+ <r (1)+ u++<7 (1)+ 

+ °12 U ( 


| n+1/2 


+ *11 «< + <r 13 ) 


C 


i) 

(2.53) 


L(i)+ U . ff (i)+ u + . (T ( 1 )+ t ,+l 

|*21 u + °22 +"23 u i|J^. t 


n+1/2 


= [< 


.(!)- 

21 


u -f a, 


(!)-.,+ 


22 U ( + *23 


(l)- u +l" 

23 ”-l(y.t;i.2) 


(c) Eq. (2.12) with £ = 3: 


(2.54) 


[<4l )+U + ^32 )+U C + *33 +U v] i k 


1 n+1/2 


= [< 


.(!)- 

31 


11 4- <T 


( 1 )” n + 4 . V 
32 ^ ■ 


(I)' 

33 


■<r 

’ J a,fc;i.3) 


1 + «( + *13 


1 n +l 

n+1/2 


1 n+l /2 


(d) Eq. (2.13) with / = 1: 

l J J »* 

= [*n )_ « 

(e) Eq. (2.13) with £ = 2: 

f <T ( 2 )+ 11 _l /r( 2 ) + ti + -4- <7® + ti + l n+ 

[*21 * + *22 U C ' *23 *17^. ^ 

r (2>- , (2)- + , (2)- +i n+1/2 

= [*21 * + *22 *C+°23 <J 0iJb;2> . 

(f) Eq. (2.13) with £ = 3: 

[-(2)+-. i 0.(2)+-,+ 4- </ 2 ^ + u + l n+1 

*31 * + *32 U ( + *33 *17 I . 

L j J ,K 


(2.55) 


i) 

(2.56) 


1 n+l 


2 ) 

(2.57) 


= tog- 


V 4- <7. 


(2) ~u + +<r' 
32 *C + * 


(2)- u +l n+1/2 

33 ” J 0,t;2.3) 


(2.58) 

Here (j, ib,n + 1/2) 6 fii is assumed in Eqs. (2.53)- 
(2.55); while (j, ib,n + 1) € is assumed in 

Eqs. (2.56)-(2.58). Also, to simplify notation, in the 
above and hereafter we adopt a convention that can 
be explained using the expression on the left side of 
Eq. (2.56) as an example, i.e., 


12 *c + * 


[<r[ 2)+ u 

= K 2 i )+ «?r + + wwn 1 ] 


.<*>+.,+ j. „<?>+*+] n+i 


According to Eqs. (2.33)-(2.35), and 

* 13 ^ contain a common factor (1 — — v^). Similarly, 

each of three consecutive pairs of coefficients defined 
in Eqs. (2.36)-(2.50) also contain a common factor. As 
a result, one concludes that: 

(a) Eq. (2.53) is satisfied by either 1 — = 0 or 

r , .T n+1/2 

[tl + (1 + l/()«^ + (1 + Vt,)u$\. k 

= [.-(i + n )«/-(i + r„)<]" Mi) 

(2.59) 

(b) Eq. (2.54) is satisfied by either 1 + vq = 0 or 

r • n+ 1/2 

[“ - (2 - *0“* + ( 1 + v i) u t j ; . k 
= [. + (S-. { )«J-(l +*)«+]" 

(2.60) 

(c) Eq. (2.55) is satisfied by either 1 + v n = 0 or 

r , ,i n+ 1/2 

« + (1 + i/()«f - (2 - v„)v + 1 

L 

= [„- (1 + ,.()«+ + (2 

(2.61) 

(d) Eq. (2.56) is satisfied by either 1 + uq + v n = 0 or 

[« - (i - 1/ ( )«^ - (i - ] ; . k 

r i n+1/2 

= [u + (1 - .()»/ + (1 - K,)«j] 

(2.62) 

(e) Eq. (2.57) is satisfied by either 1 — V( = 0 or 

r i ,in+i 

[« + (2 + - (1 - *'„)ti + t 

r i n+1/2 

= [u - (2 + i/ ( )u+ + (1 - i/„)u+] yt22) 

(2.63) 

(f) Eq. (2.58) is satisfied by either 1 — v n = 0 or 

r , , i n +i 

[u - (1 - 1/ C )t£ + (2 + ^,)< 

r i n+1/2 

= [« + (l-. ( )»+-(2+^)«f] 0iW) 

(2.64) 

Here (j, k, n + 1/2) € Qi is assumed in Eqs. (2.59)- 
(2.61); while (j,k,n + 1) € Cl 2 is assumed in 
Eqs. (2.62)-(2.64). The current 2D a scheme will 
be constructed using Eqs. (2.59)-(2.64) instead of 
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Eqs. (2.53)-(2.58). Note that Eqs. (2.59)-(2.64) im- 
ply Eqs. (2.53)-(2.58). However, the reverse is false 
unless one assumes that 


[1 - (i/ c + v v f] (1 - if) (1 - V*) # 0 (2.65) 

Note that the expressions within the brackets in 
the first three equations in Eqs. (2.59)-(2.64), respec- 
tively, can be converted to those in the last three by 
reversing the and a — ” signs. 

Let 8^\ 8^\ $ 3 ^, *i*\ and denote the 
expressions on the right sides of Eqs. (2.59)-(2.64), 
respectively. Then it can be shown that Eqs. (2.59)- 
(2.61) are equivalent to 


+ 1/2 
: ,k 



and 


t'c - » / n)4 1) + (i + » / <)4 1) + (! + 

* / n)4 1) ] 

- j ( 


(2.66) 

4°) 

(2.67) 

= 5 ( 


(2.68) 


where (j, i, n + 1/2) G fii- Also Eqs. (2.62)-(2.64) are 
equivalent to 



| [(i + + "->)4 2) + (i - "<)4 2) + (i - ^)4 2) ] 

(2.69) 

K)",* 1 = l (4 2) - 4’>) (2.70) 

and 

(Om 1 = |( 4 2) -4 2) ) (2.7i) 

where (j,k,n + 1) G ^2- 

For any (j, k, n) G let 


q(j,k,n) d = 


« \ 


v</ 


j.* 


(2.72) 


Let the 3 x 3 matrices Q [ k \ k = 1,2, and £ = 1,2,3, 
be the special cases (e = 0) of the more general 
defined in Sec. 3. Then Eqs. (2.66)-(2.68) can be ex- 
pressed as [2] 

f (j, k, n + 1/2) = Qi 1) ?((j, k ; 1, 1), n) 

+ Qi 1) W,k;l,2),n) (2.73) 

+ <4 1) ?(0\*; i,3),n) 


where (j, fc,n- 1- 1/2) G fli- Also, Eqs. (2.69)-(2.71) 
can be expressed as [2] 

q(j, k,n+ 1) = k; 2, 1), n + 1/2) 

+ <4 2) ? ( (j, *; 2, 2), n + 1/2) (2.74) 

+ Qi 2) «(0'.i;2,3),n+l/2) 

where (j, i,n + 1) G ^ 2 * The marching procedure in 
the 2D a scheme is formed by applying the marching 
steps defined by Eqs. (2.73) and (2.74) successively. 

The 2D a scheme has several nontraditional fea- 
tures. They are summarized in the following com- 
ments: 

(a) The 2D a scheme has the simplest stencil in each 
of their two marching steps, i.e., a tetrahedron in 
3D space-time with one vertex at the upper time 
level and the other three vertices at the lower time 
level. 

(b) Each of the conservation conditions Eqs. (2.53)- 
(2.58) represents a relation among the marching 
variables associated with only two neighboring 
SEs. This is a fundamental difference between the 
current method and other traditional methods. 

(c) It is shown in [2] that the 2D a scheme is neutrally 
stable if 

|i/ c | < 1.5, \i/rj\ < 1.5, and |i/ c + \ < 1.5. 

(2.75) 

As depicted in Fig. 11, the domain of stability de- 
fined by Eq. (2.73) is a hexagonal region in the 
i/{-i / n space. Moreover, it is also shown in [2] 
that Eq. (2.75) can be interpreted as the require- 
ment that the physical domain of dependence of 
Eq. (2.1) should fall within the numerical domain 
of dependence. 

(d) It is shown in [2] that the 2D a scheme has the 
following property, i.e., for any (j,t,n) G fi, 

q(j,k,n + l) ->q(j,k,n) as -+ 0 (2.76) 

if a x , a y , w , 6, and h are held constant. This 
property usually is not shared by other schemes 
that use a mesh that is staggered in time, e.g., the 
Lax scheme [15]. 

(e) The 2D a scheme is also a two-way marching 
scheme [2]. In other words, the same flux conser- 
vation relations Eqs. (2.12) and (2.13) can be used 
to construct the backward time-marching versions 
of the 2D a scheme [2]. 

This section is concluded with the following re- 
marks: 
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(a) the 2D a scheme is only a special case of the 2D 
a-p scheme described in [2]. It is a solver for the 
2D convection-diffusion equation 


du 

9t 


+a : 


du du ( d 2 u 9 2 u\ . 

di +ay + dy*)~° (2 ‘ 77) 


where a r , a y , and p (> 0) are constants. Note 
that this solver, as in the case of its ID counter- 
part, is unconditionally stable if a x = a y = 0. 

(b) It should be emphasized that, with the aid of 
Eqs. (2.19)— (2.22), (2.24), and (2.25), the 2D a 
scheme can also be expressed in terms of the 
marching variables and the coefficients tied to 
the coordinates (r,y). In other words, the co- 
ordinates (C,v) aie introduced solely for the pur- 
pose of simplifying the current development. The 
essence of the 2D a scheme , and the schemes to 
be introduced in the following sections , is not de- 
pendent on the choice of the coordinates in terms 
of which these schemes are expressed. 


3. The 2D a-e Scheme 

The 2D a scheme is neutrally stable and reversible 
in time. As explained in Sec. 1, such a scheme cannot 
be extended to become an Euler solver. As a result, it 
will be modified to become the 2D a-e scheme. Here e 
represents a special parameter that controls numerical 
dissipation. 

To proceed, note that the CEs used in Sec. 2 will 
not be used in this section. As a result, Eqs. (2.12) and 
(2.13) will no longer be assumed. Instead, the CEs to 
be used are 


*, n + 1/2) = f [c4%, k,n + 1/2)] U 
[CE£%, k, n + 1/2)] U [04%, k, n + 1/2) 


[3.1) 


where (j, fc,n + 1/2) € fli, and 


CE^(j, k,n + l) \CE[%, k,n+ 1)] U 

1 J (3.2) 

[C4 2 \j, k,n+ 1)] U [CE?\j, *, n + 1)] 


where (j, i, n-f 1) G ^ 2 - We shall assume that the total 
flux leaving the boundary of any new CE vanishes, i.e., 


and 


l h'-ds = 0, 

Js(CEi l )(j,k,n+l/2)) 


h* • ds = 0. 

JS(CE(*)(j t k,n+l)) 


(3.3) 

(3.4) 


Obviously, (i) i? 3 can be filled with the new CEs, and 
(ii) the total flux leaving the boundary of any space- 
time region that is the union of any new CEs will also 
vanish. 

Moreover, it can be shown that Eqs. (3.3) and 
(3.4), respectively, are equivalent to Eqs. (2.66) and 
(2.69). 

Proof: By subtracting the expressions on the right 
sides of Eqs. (2.53)-(2.55), respectively, from those 
on the left sides, and then multiplying the results by 
2wh/3 y we obtain the fluxes leaving CE^jjJb ,n -f 
1/2), CE&p(j, k, n + 1/2), and CE?\j, k, n + 1/2), 
respectively (see Appendix A in [2]). Because the 
flux leaving an interface from the CE on one side is 
the negative of that leaving the same interface from 
the CE on the other side, it is easy to see that the 
flux leaving CE^ 1 ^;, ib, n -f 1/2) (which is the union of 
CE^(;, ib, n-f 1/2), £ = 1,2,3) is the sum of the fluxes 
leaving CE^(j, Jb,n -f 1/2), £ = 1,2,3. Thus, the flux 
leaving CE^(j, ib, n + 1/2) can be obtained by sub- 
tracting the sum of the expressions on the right sides 
of Eqs. (2.53)-(2.55) from that on the left side, and 
then multiplying the result by 2wh/3. Furthermore, 
we have 




' 11 


+ *21 +*31 = 3, i = l,2 (3.5) 


and 


^> ± +4 i , )± +^ (t)± 


'12 


= O’ 


13 


'22 
+ *• 


32 


(*)± . 
23 


+ <4£ )=t = 0, *=1,2 


(3.6) 


With the aid of the above considerations, and the fact 
that the expressions on the left sides of Eqs. (2.53)- 
(2.55) are all evaluated at the same mesh point 
(j,k,n -f- 1/2), it becomes obvious that Eq. (3.3) is 
equivalent to the statement that is 1/3 of the 

sum on the right sides of Eqs. (2.5J>}-(2.55). By using 
Eqs. (2.33)-(2.41), we arrive at the conclusion that 
Eq. (3.3) is equivalent to Eq. (2.66). By invoking 
a similar argument involving Eqs. (2.56)-(2.58), and 
(2.42)-(2.50), we also conclude that Eq. (3.4) is equiv- 
alent to Eq. (2.69). QED. 

As a result, Eqs. (2.66) and (2.69) are shared by 
the 2D a scheme and 2D a-e scheme. In this section 
we shall describe how the other equations in the 2D a 
scheme i.e., Eqs. (2.67), (2.68), (2.70), and (2.71), can 
be modified such that the numerical diffusion of the re- 
sulting new scheme can be controlled by an adjustable 
parameter c. 

To proceed, for any (j, Jb, n -f 1/2) G fli, let 


ln+l/2 def 



(3.7) 
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where l = 1,2, 3. By their definitions, can t> e 

considered as the finite-difference approximations of u 
at ((j,fc;l,^),n -f 1/2), £ = 1,2,3. With the aid of 
Eqs. (2.26), and (2.29)-(2.31), Eq. (3.7) implies that, 
for £= 1,2, 3, 

%"?M) = [ u ~ 2 { u < u t + ( 3 ‘ 8 ) 

In Fig. 12(a), P, Q, and R are three points in the 
C-77-u space. Using the coordinates given in the same 
figure and Eqs. (2.51a)-(2.51c), it can be shown that 
these points are on a plane represented by 

« = K); > r /2 «-i^)+K)-t l/2 ^- 1*1)+ v'jT' 2 

(3.9) 

where 

■4r /! = 5 (««r.s + c‘.x + -ssss) «*m> 

(“()?, t 1,J = («0,w£ - ««“«?)) Mt (311) 

and 

K C ,,S = M-l (3.12) 

Equation (3.9) implies that point O depicted in 
Fig. 12(a) is also a point on the plane that contains 
P, Q , and P. Moreover, for every point on the plane 
represented by Eq. (3.9), including point O, we have 


(!) n =«>;r. - (£) -«« 


+1/2 

(3-13) 

As a result of the above considerations, 

(ti^)”^ 1 ^ 2 , and (t^)^** 1 ^ 2 can be considered as the 
finite-difference approximations of u , dti/dC, and 
3tx/9r/ at the mesh point (j, t,n -f 1/2), respectively. 
Note that Uy^ +1 ^ 2 generally is different from 2 
which is defined by Eq. (2.66). Because Eq. (2.66) is 
equivalent to the conservation condition Eq. (3.3) and 
thus a fundamental part of the 2D a-e scheme, ti*. * +1 ^ 2 
will not be used in the future development. 

To proceed, let 


(V+)"+ 1 /2 d S? \»>+l/2 

(“( >),k ~ 6 *Vi,* 

/ /+\n+l/2 def £^7/ / \n+l/2 

~ g IVj.i 

>_ 4 „) 


(3.14a) 

(3.146) 

(3.15a) 


and 

«)m 1/! = 5(«S 1> -4 1> ) (3.15*) 

where s^, and 8 3^ are the expressions on the 

right sides of Eqs. (2.59)-{2.61), respectively. Sev- 
eral comments can be made about Eqs. (3.14a,b) and 
(3.15a, b): 

(a) Eqs. (3.15a, b) simply say that the expressions on 
the right sides are denoted by the notations on 
the left sides, respectively. Contrarily, Eqs. (2.67) 
and (2.68) say that the marching variables on the 
left sides are assigned the values of the expressions 
on the right sides, respectively. 

(b) Because the expressions on the right sides of 
Eqs. (3.15a, b) are functions of the marching vari- 
ables at the nth time level, so are (ti£ + )"j^ 2 and 

(x i°+ > i f>+1 / 2 
\ U rj )j,k 

(c) Eqs. (2.30), (2.67), (3.13), (3.14a), and (3.15a) 
imply that both (u^)?* 1 ^ 2 and (u£ + )*?J^ 2 are 
numerical analogues ot du/dC at the mesh point 
(j, k, n -f 1/2), normalized by a£/ 6. Similar 
interpretations can be given to (uj*')”* 1 ^ 2 and 

K + );,t 1/2 - 

(d) (tij + )?J 1 ^ 2 and (uj*)?* 1 ^ 2 are defined as a re- 

sult of a geometric construction involving finite- 
difference approximations. Generally, numeri- 
cal dissipation will be introduced as a result of 
using such numerical analogues. On the other 
hand, (t^*)”,* 1 ^ 2 and (n® + )”^ 1 ^ 2 emerge from 
the development of the 2D a scheme which is free 
from numerical dissipation. In the following, the 
2D a-e scheme will be constructed such that (i) 
(«+ )£* is a weighted average of (u^ )JJ^ J and 
K+);* 1 ' 2 with the weight factors 2e and 1 — 2c, 
respectively; and (ii) (u+)? is a weighted aver- 
age of with the weight 

factors 2e and 1 — 2e, respectively. As a result, nu- 
merical dissipation may be controlled by varying 
the value of e. 

To proceed, for any (j, Jb,n 4- 1/2) £ fix, let 

( d«+)”t i/2 * 2 [« + )j,: i/2 - («r)?.t i/2 ] (3.16) 

and 

wit 11 * = 2 [k+) 5 ! 1/2 - «+)-r /2 ] (3- 17 ) 

Then, with the aid of Eqs. (3.8), (3.11), (3.12), 
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(3.14a,b), and (3.15a, b), one has 

(*‘<ff" = j [(«+<«? -2<)’ w) 


(3.18) 


and 




(3.19) 


Thus both (dujf)"* 1 ^ 2 and (du*)?* 1 ^ 2 are functions 
of the marching variables at the nth time level . As 
will be shown shortly, they play a key role in the first 
marching step of the 2D a-e scheme. 

Next we consider Fig. 12(b). For any (;, Jfc, n-f 1) € 
Q 2 » let 

/ a* \ n+l/2 

U G\*;2,0 = ( u + T Ut ) ( 3 - 2 °) 

where £ = 1, 2,3. Then it can be shown that 

1 n+l/2 


(Off = ^«>ff (3.264) 

(«; + )ff = j(4’ > -*‘ ,) ) ( 3 - 27 «) 

and 

K + )".* 1 d " | (4 2) - * i 2> ) (3-276) 

where s^\ s^\ and are the expressions on the 
right sides of Eqs. (2.62)-(2.64), respectively. Because 
these expressions are functions of the marching vari- 
ables at the (n -f l/2)ti time ievei, so are 

“<* K + )i,t 2 * 

With the above preparations, the current coun- 
terparts to Eqs. (3.16) and (3.17) sire 

(<*«<)]? = 2 - K +)# 1 ] (3.28) 


and 


('Off = 2 [(Off - (Off ] o-m) 


respectively, where (j,k,n + 1) € & 2 - With the aid 
of Eqs. (3.21) and (3.24)-(3.27a,b), Eqs. (3.28) and 
(3.29) imply that 


. , , r / \ in+i/2 (dut)^ 1 = i[(u + 2ut + 2u+) 


In Fig. 12(b), the points P, Q, and R are on a 
plane represented by 

“ = (Off « - 30+ (“iff (i - *“i) + «ff 1 

(3.22) 

where 

2. / ti /n +l I _ t /n+l i , # /T »+l ^ /o oo\ 

i.* ~ 3 V*(i,*;2,i) + tt 0\*;2,2) + u U,k,2,3)) {o.Z6) 


(Off - («'ffw) - <©«>) /+<• ( 3 - 24 > 


- («-k+ 2 <)"*XJ 


(3.30) 


and 


W = 5 [(* ■ + 2n ? + 2 <CX> 

-(« + 2.+ -4 <)^"J 


(3.31) 


and 


— ( u J!it;2,3) “ u 0\M,i)) / AT? * (3.25) 

Moreover, the current counterparts to Eqs. (3.14a,b) 
and (3.15a,b) are 

(Off = T«ff ( 3 - 26 «) 


Thus both (du^V?J 1 and (dti*)? J 1 are functions of the 
marching variables at the (n -f 1/2) th time level . 

The 2D a-e scheme can now be stated using the 
above definitions. It consists of two marching steps. 
The first is formed by Eq. (2.66), 

(Off" = (Off" + <('Off " < 3 - 32 ) 


and 


Kff' ! =H"ff'’+‘K)i; 


+ \n+l/2 


(3.33) 


where (j, k, n+ 1/2) € fir- It was explained earlier 
that the expressions on the right sides of Eqs. (3.32) 
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and (3.33) are functions of the marching variables at 
the nth time level. Moreover, according to Eqs. (3.16) 
and (3.17), Eqs. (3.32) and (3.33) can also be expressed 
as 


K + )”t 1/2 = K + )J i t 1/2 +(e- 1 /2)(dt,+);+ 1/2 (3.34) 


and 


«)-: 1/2 = K + ff /J + ( £ " 1 / 2 )K)5 1/2 (3-35) 


respectively. 

The second marching step is formed by Eq. (2.69), 


u 


+ \n + 1 

c 


= K + )zr 




(3.36) 


and 


M)S‘ = (OS' +«)”.{' 


(3.37) 


where (j, Jb,n -f 1) € K was explained earlier that 
the expressions on the right sides of Eqs. (3.36) and 

(3.37) are functions of the marching variables at the 
(n -f l/2)th time level. Furthermore, according to 
Eqs. (3.28) and (3.29), Eqs. (3.36) and (3.37) can also 
be expressed as 

K);,r = « + )-,r + (* - i/ax**?)# 1 . (3.38) 

and 

«)",{' = + (« - (3.39) 

respectively. 

At this juncture, note that: 

(a) With the aid of Eqs. (3.15a,b) and (3.27a, b), it is 
seen that Eqs. (3.32), (3.33), (3.36), and (3.37), 
respectively, are reduced to Eqs. (2.67), (2.68), 
(2.70), and (2.71) when e = 0. As a result, the 
2D a-e scheme becomes the 2D a scheme when 
e = 0. 

(b) For the special case with e = 1/2, Eqs. (3.34), 
(3.35), (3.38), and (3.39) are reduced to the 
forms that represent the finite-difference approx- 
imations defined in Eqs. (3.11), (3.12), (3.14a,b), 
and (3.24)-(3.26a,b). However, Eqs. (2.66) and 
(2.69), which are independent of e and therefore 
always part of the a-e scheme, are the results of 
the flux conservation conditions Eqs. (3.3) and 
(3.4). 


(c) With the aid of Eqs. (2.30) and (3.18), Eq. (3.32) 
can be rewritten as 


(u ,) n+1/2 - _f«°+'l n+1 ' 2 
- aC ^C t),k 

+ 5[(S +4 “'-^)’ tt . w) 


(3.40) 


^ (0 ("fW *** K)o\*;i,2) be 

identified with the values of ti, du/d£ and du/drj 
at the mesh point ((j, k ; 1, 2), n), respectively; and 

( fi ) («c)&.m.i) and be iden ‘ 

tified with the values of ti, du/dC and du/dr) 
at the mesh point ((j, it; 1, l),n), respectively. 
Then it can be shown that the expression within 
the brackets on the right side of Eq. (3.40) is 
0 (aC,at 7). Furthermore, because Eq. (2.28) is 
applicable within SE(j,Jb,n) only, the expression 
that is enclosed within the first bracket on the 
right side of Eq. (2.28) is 0 (a£,a<). From the 
above considerations, one concludes that the ad- 
dition of the extra term involving e on the right 
side of Eq. (3.40) may result in errors that are sec- 
ond order in a£, atj, and a f. In other words, the 
addition of the term involving e does not result 
in a scheme of lower order of accuracy. A similar 
conclusion is also applicable to Eqs. (3.33), (3.36), 
and (3.37). 


(d) Equations (3.11), (3.13) and (3.14a) imply 

that (ti^)?^ 1 ^ 2 is proportional to the directional 
derivative along the C-direction on the plane that 
contains points P, Q, and R which are depicted in 
Fig. 12(a). According to Eq. (3.16), (du^)J^ 1/r2 

is twice the difference between (ti£ + )£^ 2 and its 
counterpart in the 2D a scheme. Note that the 
variable (dti x )£, that appears in Eqs. (3.2) and 
(3.10) of [1,9], plays a role in the ID a-e scheme 
[1,9] similar to that of (du^V^J 1 ^ 2 in the present 
2D a-e scheme. It can be shown that (du x )” is 
equal to the difference between two slopes. The 
first slope is the central difference given on the 
right side of Eq. (3.10) in [1,9]. The second slope 
is the counterpart of the first in the ID a scheme. 
Thus the 2D a-e scheme is a natural extension of 
the ID a-e scheme. 

The a-e scheme will take the form of Eqs. (2.73) 
and (2.74) if the 3x3 matrices Qi k \ k = 1,2, and 
£ = 1,2,3 take the form to be specified immedi- 
ately. Because of the limitation imposed by the cur- 
rent double-column format, their elements can only 
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be specified row by row, sequentially. Furthermore, 
for the sake of simplicity, the elements of the matrices 
ZQ^ will be listed: 

3Q ( 1 1) : 1 - v ( - v -(1 - - i/„)(l + v ( ), -(1 - vc - 

i/,)(l + *>,), 1 - e, -(1 + - 2e), -(1 + u„- 2e), 

1 — e, — (1 + — 2c), — (1 + Vr) — 2c). 

3Q ( 2 1) : 1+Hf, (1 +j/()(2-i/ c ), -(l+V()(l+i^), -(1-0. 
—(2 — - 4e), 1 + v n - 2c, 0, 0, 0. 

3 Q^: 1 + Vff t -(1 + i/„)(l + i^c), (1 + v v ){2 - v n ), 0, 0, 
0, -(1 - e), 1 + 1/ ( - 2e, -(2 -u„- 4c). 

3Q ( x 2) : 1 +i/ ( +i/,,(1+i/ ( +i/,)(1-i' C ), (l+*'c+i',,)(l-i'i,), 
-(1 - e), -(1 - 1 / ( - 2c), -(1 -Vr,- 2c), -(1 - e), 
— (1 — •'C — 2c), —(1 — v n — 2c). 

3Q ( 2 2) : 1 - »/ ( , -(1 - * c )(2 + */ ( ), (1 - ^)(1 - i/,), 1 - c, 
-(2 + i/(- 4c), 1 - u v - 2c, 0, 0, 0. 

3Q ( 3 2) : 1 - i/„, (1 - ^)(1 - v ( ), -(1 - I4j)(2 + i/q), 0, 0, 
0, 1 - c, 1 - f/ c - 2c, -(2 + ^ - 4c). 

It is shown in [2] that (i) the 2D a-c scheme is 
unstable ifc<0orc> 1, and (ii) numerical diffusion 
increases as c increases, at least in the range of 0 < 
c < 0.65. In order to suppress numerical oscillations 
near a discontinuity, one may be forced to choose a 
large c. However, with such a choice, the smooth part 
of a solution may become highly diffusive. To resolve 
this dilemma, in the following, we shall construct a 
generalization of the a-c scheme. 

To proceed, let (j, fc,n + 1/2) E fii and con- 
sider Fig. 13(a). This figure is essentially identical to 
Fig. 12(a) except that point O in the latter is replaced 
by point O* in the former. The coordinates of point 
O* are (jaC, k&rj, u”^ 1 ^ 2 ) where is defined in 

Eq. (2.66). Thus point O * generally is not on the plane 
that contains points P, Q, and R. Let planes #1, #2, 
and #3, respectively, be the planes containing the fol- 
lowing trios of points: (i) points O*, Q, and R ; (ii) 
points O *, P, and P; and (iii) points O * , P, and Q. 
Then in general these planes differ from one another 
and from the plane that contains points P, Q , and R. 
In the following, first we shall study the former three 
planes. 


(«S 2) )"£ 1/2 = (2*i + *3 )/a< (3.44) 

(“£°)?,t 1/2 = f (*i-* 3 )M»? (3-45) 

(“S 3) 8 ,/5 = (*i"*j)K (3-46) 

(li ( 3) )n+l /2 <& (2xi + X2 y Aq ( 3 .47) 

Moreover, for £ = 1,2,3, let 
(„,, ,)»♦'« «(„(<>)»«/> §+(„(<>>;«'’ g (3.48) 


and 


and 


(4°)7i 1/a = («c°);.£ 1/2 f + « ) >;,£ 1/2 ( 3 -49) 

With the aid of Eqs. (2.20) and (2.22), we have 

(4Xt‘ p = (“(‘’l/.t 1 ' 1 + 3 ; 1,1 (»•»>) 

and 

/.//Vn+l/2 («> + 6)AC ,(Osn+l/2 


2wh 

(w - 6)at? / (Qx»+ 1/* 
2u)A 1 " <#•* 


(3.51) 


Combining Eqs. (3.42)-(3.47) with Eqs. (3.50) and 
(3.51), one has 


(4>)) ? « , . = -±(x, + l3 ) 


, (i)\ n +i/2 _ (35 4- W)X2 + (35 ~ 1*0*3 

'i> k 2wh 




(3.52) 

(3.53) 

(3.54) 


a preliminary, let 

(3.41) 

, f 2 l\n+l /2 (36 + tl/)*l + 2«;*3 

^ '>•* 2tt)/» 

(3.55) 

(„(.>).«« ^_ (2IJ+I3)/ac 

(3.42) 

r«( 3 )i?t 1/2 - ^£1 

y e } }> k ~ 2w 

and 

(3.56) 


(3.43) 

, f3K«+l/2 _ (» ~ 35)ll + 2tl>* 2 
' u v 'jy* 2 tvh 

(3.57) 


15 


With the above preparations, it can be shown 
that, for each £ = 1,2,3, plane #£ is represented by 

« = 1/2 (C - JAC) + «))£ l/2 (r, - kAtj) 

. n+l/2 

+ u j,k . 

(3.58) 

if the coordinates (C,»?) are used. Alternatively, it is 
represented by 

» = 1,2 (* - *M> + (r - «,*) 

(3.59) 

if the coordinates (x,y) are used. Using Eqs. (3.58) 
and (3.59), one concludes that, at any point on plane 
#£, £ = 1,2,3, we have 


«.d (£) 


+1/2 

k 


and 


(3.60) 




^du 

(3.61) 

Note that Eq. (3.60) is the current counterpart of 
Eq. (3.13) which is applicable to any point on the plane 
that contains points P> Q, and R. Let Vti be the gra- 
dient of tx. Then Eq. (3.61) implies that, at any point 
on plane #£, £ = 1,2,3, we have 


iv«i=(*/);,£ 1/2 = f 


y(«? ) ) 2 +(4 <) ) 2 


n+l/2 




(3.62) 

To proceed further, we introduce the current 
counterpart of Eqs. (3.14a,b), i.e., 

(„W+)?+ 1/2 « ^(t£>)"£ 1/2 (3.63a) 


and 


<*•«») 

Then E<|s. (3.11), (3.12), (3.14a, b), and (3.41)-(3.47) 
imply that 

(u' + )7t 1/2 = l [« ( ( 1)+ + «< 2)+ + «4 3)+ £ 1/2 (3.64) 
and 

(os 1 " = | k ,h + + » r]"r (»«) 


i.e., (i) is the simple average of tc^ + , £ = 1,2,3. 

and (ii) uj+ is the simple average of ui^ + , £ = 1,2,3. 
Note that, for simplicity, hereafter we may suppress 
the space-time mesh indices if no confusion could oc- 
cur. 

The first marching step of the generalized a-e 
scheme will be formed using Eqs. (2.66), (3.34), and 
(3.35) except that (i) in Eq. (3.34) is replaced 

by a nonlinear weighted average of ti^ + , £ = 1,2,3, 
and (ii) in Eq. (3.35) is replaced by a nonlinear 

weighted average of £ = 1,2,3. The design of 

these weighted averages will be guided by the require- 
ment that the weight assigned to a quantity associated 
with plane # £ is greater if 0* is smaller. This require- 
ment is similar to that used in the construction of an 
ID Euler solver described in [1,9]. 

To proceed, for any a > 0, we shall define u^ + 
and , the weighted-average counterparts of u'+ and 
u'+ f respectively. Let 


u“ + = 0, if 0! = 0 2 = 03 = 0; (3.66a) 


and 


w+ (M 3 ) a t4 1H + (Mo a t4 2)+ + ( wtff * 

“ c (M 2 ) a + (M 3 )° + (Mi)° 

otherwise. Also, let 


(3.666) 


< + = 0, if = 0 2 = 03 = 0; (3.67a) 


and 

,+ (020 3 ) a 4 1)+ + (Mi) a tff )+ + ( w«4 3H 

(M 2 )° + (0 2 0 3 )° + (Mi)“ 

(3.676) 

otherwise. Note that the denominators of the fractions 
on the right sides of Eqs. (3.66b) and (3.67b) vanish 
if or > 0, and any two of 0\, 02, and 0 3 vanish. Thus, 
consistency of the above definitions requires the proof 
of the proposition: 0\ = 0 2 = ^3 = 0, if any two of 0 
0 2 , and 03 vanish. 

Proof: As an example, let 0 X = 0 2 = 0. Then 
Eq. (3.62) implies that u * ^ = Uy 1 ^ = 0, £ = 1,2. In 
turn, Eqs. (3.52)-(3.55) imply that x\ = x 3 = *s = 0. 
0 3 = 0 now follows from Eqs. (3.56), (3.57), and (3.62). 
QED. 

As a result of Eq. (3.66b), we have 


= 


“C ’ 
,( 2 )+ 
“c ’ 

,(3)+ 


if 0i = 0, 0 2 > 0, and 0 3 > 0; 

if 0 2 = 0, 0i > 0, and 0 3 > 0; 

if 0 3 = 0, 0i > 0, and 0 2 > 0. 

(3.68) 
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Assuming Ot > 0, £ = 1,2,3, we have 


and 


+ wwyr + (w a «4 : 


(2)+ 


,.( 3 )+ 


U, 1 = 


(i/0i) a + {i/e 2 ) a + (i MO" 


(3.69) 


Thus the weight assigned to u^ + is proportional to 
( l/9 t ) a . By using Eqs. (3.64), (3.66a), and (3.69), one 
arrives at the conclusion that 


«”+ = «'+, if 8 l = 6 2 = 0 3 . (3.70) 

Obviously Eqs. (3.68)-(3.70) are still valid if each sym- 
bol C is replaced by the symbol rj. 

On the smooth part of a solution, 0i, 0 2 , and 63 
are nearly equal. Thus the weighted averages tx™ + and 
are nearly equal to the simple averages and 
tt'+, respectively (see Eqs. (3.64) and (3.65)). In other 
words, the effect of weighted-averaging generally is not 
discernible on the smooth part of a solution. 

Next let 1) £ Q 2 and consider Fig. 13(b). 

The third coordinate tx?* 1 of point O* is that defined 
in Eq. (2.69). Let planes #1, #2, and #3, respectively, 
be the planes containing the following trios of points: 
(i) points O*, Q, and R\ (ii) points O * , P, and P; and 
(iii) points 0*, P, and Q. In the following, we shall 
study these planes. 

As a preliminary, let 


Vt 


— f tx ,n+1 — tx n+1 / — 1 9 3 


(3.71) 


(t|( 1) )»+ 1 H f (2y 2 + y3 )/AC 

(3.72) 

« ) )i,t 1 = (y2 + 2y 3 )/A I , 

(3.73) 

(< 2) )",t ld = -(2y! + y3)/AC 

(3.74) 


K 2) )"t ld = (ya-yO/A'J (3.75) 


(«( 3) )",t ld = (y 2 -yi)/A< (3.76) 

and 

(^ 3) }”.t 1 d = ~(2yi + y 2 )/A»7 (3.77) 

Moreover, for £ = 1,2, 3, let 




(3.78) 


( W = gj + (« W S (3 79) 




With the above definitions, Eqs. (3.50) and (3.51) re- 
main valid if each upper index n -f 1/2 is replaced by 
n + 1. As a result, Eqs. (3.72)-(3.77) imply that 


(«* ) )”,t 1 = 2^(l/2 + Ife) (3.80) 


(««>”» 1 = . 

(36 + w)y 2 + (36 - w)y 3 

(3.81) 

2 wh 


w— £ 

(3.82) 


_ (36 + U))yi + 2 uij/ 3 
2u/h 

(3.83) 


( u (3))n+l _ _3yi 
V x )j,k 2t0 

(3.84) 

and 



< W • 

(to - 36)yi + 2wt/2 
2 wh 

(3.85) 


With the above preparations, the earlier develop- 
ments that involve Eqs. (3.58)-(3.70) can be repeated 
for the current case with the only change being the re- 
placement of each upper index n-fl/2 by n+1. Partic- 
ularly, («c + )?,t 1 ( u v + )jjs l can be defined using 
Eqs. (3.66a, b) and (3.67a, b) with the understanding 
that each symbol in these equations is associated with 
the mesh point (j, Jb,n-b 1). 

The generalized a-e scheme, referred to as the 
weighted-average a-c scheme, can now be stated. It 
consists of two marching steps. The first is formed by 
Eq. (2.66), 

(«< + )m 1/2 = (u»+);+ 1/2 +( C -l/2)(du+); i t 1/2 (3.86) 
and 

«)”t 1/2 = (u“+);; 1/2 +(e-l/2)(du+)"+ 1/2 (3.87) 

where (j,k, n + 1/2) G fi). The second is formed by 
Eq. (2.69), 

(«?&' = (“Djt 1 + (<■ - maw? < 388 > 
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and 

«) "t 1 = K + )"t 1 + (e - l/2)(d«+)?+ 1 (3.89) 

where (;,fc,n + 1) € 0 2 . 

Note that, according to Eq. (3.62), evaluation of 
(i 9i) Q does not involve a fractional power if a is an 
even integer. Because a fractional power is costly to 
evaluate, use of the generalized a-e scheme is less costly 
when a is an even integer. 

4. The 2D Euler Solvers 


are the space-time mass, i-momemtum component, y- 
momemtum component, and energy current density 
vectors, respectively. 

As a preliminary, let 

r m ,t = drjd u t , and fa/=dfa/d« t (4.13) 

where = 1,2, 3,4. The matrices formed by f^ t 

and m,e = 1,2, 3, 4, respectively, are given in [2]. 

Because /* and /£, m = 1,2, 3, 4, are homogeneous 
functions of degree 1 [19] in tii, ti2, tf3> and U4, we have 


We consider a dimensionless form of the 2D un- 
steady Euler equations of a perfect gas. Let />, ti, v, p, 
and 7 be the mass density, z-velocity component, y- 
velocity component, static pressure, and constant spe- 
cific heat ratio, respectively. Let 

ui = p, u 2 = pit, 1 i 3 = pv (4.1a) 

«4 = p/(7 - 1) + p(« 2 + n 2 )/2 (4.16) 

/? = « 2 , 

/f = (7 ~ 1)«4 + (3 - 7)(ti 2 ) 2 /(2u!) 
-(7-l)(« 3 ) 2 /(2«i) 

/| = «2 « 3 /«l. 
fl = 7«2«4/«l 

- (!/2)(7 ~ 1)«2 [(« 2 ) 2 + ( U3 ) 2 ] /(« i ) 

/? = «3, 

fl 

/I = (7 - 1)«4 + (3 - 7)(« 3 ) 2 /(2u!) 

-( 7-l)(ti2) 2 /(2m) 

and 

/» =7u 3 u 4 /« 1 

- (l/2)(7 - 1)«3 [(«2) 2 + («3) 2 ] /(«l) 2 
Then the Euler equations can be expressed as 

dfa d fa _ 

-dr + -te + -di~ 0 ' m - 1 ' 2 ’ 3 ’ 4 


(4.2) 

(4.3) 

(4.4) 


2 (4-5) 


(4.6) 

(4.7) 

(4.8) 


/ ' 
JS(V) 


h m ■ ds = 0, m = 1,2, 3, 4 


where 


fa = J2 fm,i «/. “d fa = Yj ut ■ ( 4 - 14 ) 


1=1 


1=1 


For any (x,y,t) 6 SE (j,k,n), v m (x,y,t), 

fa(x,y,t), fa(x, y,t), and h m (x,y,t), respectively, 
are approximated by tz£,(a:,j 

fa(x,y,t;j,k,n), and 

h^ n (x, y, t ;j, k, n). They will be defined shortly. For 
m = 1,2, 3, 4, let 


ti m(x,y,t;j,k,n) d = (u m )",t + («mr)",i(* - *;,i) 


+ (« m y)" it (y - ! Ij,ic) + (wm«)",t(f - f n ) 


(4.15) 


where (v m )? t , (tw)" it , («my)" it . and are 

constants in SE(j, Jb,n). Obviously, they can be con- 
sidered as the numerical analogues of the values of u m , 
dum/dx , dum/dy , and dum/dt at (xj,y k ,t n ), respec- 
tively. 

Let (fa)l k , (fa)lt, {fa,i)lu, and (fa^ de- 
note the values of fa, fa, fa t , and fa t , respectively, 
when Mfn ) m — 1,2, 3, 4, respectively, assume the val- 


Am = (/*,/*, « m ), m = 1,2, 3,4 (4.12) 



ues of (u m )" *, m = 1,2, 3,4. Let 


(4.9) 

us ,,)?.* d =ix 
/= 1 

(4.16) 



(4.17) 

(4.10) 

/=1 

4 


’ 3 is 


(4.18) 

Z=1 

4 


(4.11) 

( A * " DOW ”*)".* 

(4.19) 


*=1 

4 


(4.12) 

/=1 

(4.20) 
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and 


VDj.k = ( 4 - 21 ) 

Z=1 


Because (i) 


df X m 

dx 


_ v ' f * 

- 2-, o_ » 


*=1 


m = 1,2,3,4 (4.22) 


and (ii) the expression on the right side of Eq. (4.16) 
is the numerical analogue of that on the right side 
of Eq. (4.22) at (x/,y*,* n ), (/m*)”,fc can be consid- 
ered as the numerical analogue of the value of df^/dx 
at ( Xj,y k ,t n ). Similarly, (/^ y )" >t , (/mr)",t» 

( fmt)j,k can considered as the nu- 
merical analogues of the values of dfa/dy, df^/dt, 
dfm/dx, df^/dy, and df^/dt at(x,, y k ,t n ), respec- 
tively. As a result, we assume that, for m = 1,2, 3, 4, 

/£(*,»,* ;i. *,«) = r (£)£* + (/L)j, t (* - *i.») 




(4.23) 


and 


/£*(*, y,t ;j,Jb,n) d ^ f (/»)£* + («*(* - 


+ (^y)?, i (y-yi,t) + (« i (^-n 


(4.24) 


Also, as an analogue to Eq. (4.12), we assume that 
^(x,y,t;j,k,n) = f (/£*(*, y,i,j,k,n), 


fm (*> y. < ; j, *> «), <»(*> y. < ;j> M)) 


(4.25) 


Note that, by their definitions: (i) /£, /■£,, fm,t> and 
fmi’ = 1> 2 ,3, 4, are functions of («*»)?,*, m = 
1,2, 3, 4; (ii) (/* r ); >t and (/»«)£*, m = 1,2, 3,4, are 
functions of («„»)”* and (ti m r)" |t , m = 1, 2,3,4; (iii) 
(fmy)j,k (fmy)j,k’ ni = 1,2* 3, 4, are functions of 
(u m )? t and (umy)" it , m = 1,2, 3, 4; and (iv) (/£,)",* 
and (fmt)j,k are functions of («m)” it and (u mt )? k , m = 
1,2, 3, 4. 

Moreover, we assume that, for any ( x,y,t ) € 
SE(i, k, n), and m = 1, 2, 3, 4, 


au^(z,y,t;j,fc,n) g/”(»,y,f ;j,fc,n) 
dt dx 


(4.26) 


Note that Eq. (4.26) is the numerical analogue of 
Eq. (4.10). With the aid of Eqs. (4.15), (4.16), (4.20), 


(4.23), and (4.24), Eq. (4.26) implies that, for m = 

1 , 2 , 3 , 4 , 




1=1 3 ’ 


(4.27) 


Thus (««,<)",* are functions of (u m )£ t , (u mx )" (t , and 
(u m y)”f From this result and the facts stated fol- 
lowing Eq. (4.25), one concludes that the only inde- 
pendent discrete variables needed to be solved in the 
current marching scheme are (tim)" ik , (tw)”,*, and 
( u my)j t k- 

Consider the conservation elements depicted in 
Figs. 7(a) and 8(a). The Euler counterparts to 
Eqs. (2.12) and (2.13), respectively, are (i) 

/ h' m ds = 0 (4.28) 

J S(C Eg (j ,k ,n+l/2)) 

where (j, Jb,n -f 1/2) € fii; and (ii) 


f (3) 0 


(4.29) 


where (j, Jb,n-f 1) € ^ 2 - Note that £ = 1,2,3, and 
m = 1,2, 3, 4 in Eqs. (4.28) and (4.29). 

Next we shall introduce the Euler counterparts 
of Eqs. (2.24), (2.25), (2.30), and (2.31). For any 
O', *, ») € ft, let 



where m,£ = 1,2, 3, 4. The normadized counterparts 
of those parameters defined in Eqs. (4.30) and (4.31) 
are 


«>?> = 

(4.32a) 


(4.326) 

d — (f< 

~~ 2a£ 

(4.33a) 

^ \n 

(4.336) 
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In the following development, for simplicity, we 
may strip from every variable in an equation its in- 
dices i, Jb, and n if all variables are associated with 
the same mesh point (j,fc,n) £ Q. Let F C+ and F T?+ , 
respectively, denote the 4 x 4 matrices formed by 
and m,£ = 1,2, 3, 4. Let I be the 4 x 4 identity 
matrix. Then the current counterparts to Eqs. (2.33)- 
(2.50) are 


Eg* d = I - F <+ - F> + 

(4.34) 

Eg* = f ±(7 - F c+ - F V+ )(I + i* + ) 

(4.35) 

Eg* d = ±(7 - F (+ - F r,+ ){I + F” + ) 

(4.36) 

Eg*H f 7 + F< + 

(4.37) 

Eg* d = T(7 + F c+ )(27 - F <+ ) 

(4.38) 

Eg* = f ±(7 + F c+ )(7 + F n+ ) 

(4.39) 

Eg* = f 7 + F v+ 

(4.40) 


and 3 can be converted to their counterparts in this 
section. The latter will be referred to as the Euler 
images of the former. 

Equations (4.28) and (4.29) are evaluated in Ap- 
pendix B of [2], Let u, u + , and £f+, respectively, be 
the 4x1 column matrices formed by ti+£, and 
u mrjy m = 1,2, 3, 4. Then the results can be expressed 
as: (i) 


i n+l/2 


= [Eg>-«+ tig-t ? + w/) < 4 - 52 > 


where (j, k,n + 1/2) € Di and l = 1, 2, 3, and (ii) 

[Eg»*+isr< t+^Z 1 

- fE (2)- « + E (2)- t?+ + E (2)_ u + l n+1/2 ( 453 ) 


Eg* = f ±(i + r ,+ )(/ + f (+ ) 
Eg* = f T(I + F'+XM - F» + ) 
Eg* <y.f j + F (+ + p^+ 
Eg* T (J + F (+ + F" + )(7 - F< + ) 
E&* d = ^(7 + 7*+ + F* + )(7 - F” + ) 
Eg* d = J - F c+ 

Eg* d = ±(7 - F <+ )(27 + F <+ ) 
Eg* = f - F c+ )(7 “ F* + ) 

^2)i j __ + 

Eg* = f =f(7 - F V+ )(I - F <+ ) 


(4.41) 

(4.42) 

(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 


where (j, k,n + 1) E ^2 and £ — 1,2,3. Note that the 
six equations given in Eqs. (4.52) and (4.53) are the 
Euler images of Eqs. (2.53)-(2.58), respectively, under 
the substitution defined by Rule 2 given above and 
Rule 3: u, u t , t£, and u+ be replaced by u, u«, tT ( + , and 
u + , respectively. 

Note that an Euler image is a matrix equation. 
Because (i) matrix multiplication is noncommutative, 
and (ii) the coefficient matrices of an Euler image are 
mesh-point dependent, an Euler image is more difficult 
to deal with than its counterpart in the 2D a or a-c 


(4.48) scheme. 

Because the Euler images of Eqs. (3.5) and (3.6) 
(4*49) ajg also true [2], by summing over the three equa- 
tions given in Eq. (4.52), one concludes that, for any 
(4 ' 50) (j, k,n + 1/2) € Oi, 


and 


Eg* = ±(7 - F” + )(27 + F* + ) (4.51) 

Note that Eqs. (2.33)-(2.50) become Eqs. (4.34)- 
(4.51), Eqs. (4.34)-(4.51), respectively, under the fol- 
lowing rules of substitution: 

Rule 1: 1, i/(, and i be replaced by /, F^ + , and F* 7 *, 
respectively 

Rule 2: and <rg* be replaced by Egfe and Egfe, 

m,£ = 1,2,3, respectively. 

As will be shown, under the above and other rules of 
substitution, many other equations given in Secs. 2 


«sr*- 


1 V feW - ,7 4- E (1)_ u + 4- E (1)_ u + 1 ” 

j 2^ fen u + ‘2 «( + J 0 . (fc;M) 


(4.54) 


As a result, u **jfe 2 can be evaluated in terms of the 
marching variables at the nth time level. Similarly, by 
using Eq. (4.53), one concludes that, for any ( j,k,n + 
1) G 0 2 , 


»7T = 




1=1 


n+ 1/2 (4.55) 
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As a result, tP^ 1 can be evaluated in terms of the 
marching variables at the (n 4- l/2)th time level. 

For any (j, Jfc,n -f 1/2) £ fii, the matrices 

( E ml + )" l t 1/2 . m ’ £ = 1,2,3, are functions of t£| 1/2 . 
Thus they are also functions of the marching variables 
at the nth time level. Assuming the existence of the in* 
verse of each of the matrices (E^ + )?* ^ 2 , m = 1,2, 3, 
one can define 




(4.56) 


where £ = 1, 2,3, and the inverse of a matrix A is de- 
noted by [A]" 1 . As a result of their definitions, s£ l \ 
£ = 1,2,3, can be evaluated using the marching vari- 
ables at the nth time level. 

For any (j, fc, n+1) £ ^ 2 , the matrices (E^ + )”£ 1 , 
m, £ = 1, 2, 3, are functions of fi? J 1 . Thus they are also 
functions of the marching variables at the (n -f l/2)th 
time level. Assuming the existence of the inverse of 
each of the matrices (E^*)?* 1 , m = 1,2,3, one can 
define 


x [Eg>-«+ X&-S+ + 


(4.57) 


As a result of their definitions, S} 2 \ £ = 1,2,3, can 
be evaluated using the marching variables at the (n -f 
l/2)th time level. 

With the aid of Eqs. (4.34)-(4.51), (4.56), and 

(4.57), Eqs. (4.52) and (4.53) imply that [2] 

[« + (/+ F<+) u+ + (7 + F v+ ) t?+]J^ /2 = 5 X (1) , 

(4.58) 

[u - (27 - i* + ) «+ + (7 + F” + ) ^ + ]”^ /2 = 5 2 (1) , 

(4.59) 

[«+ (7+ F< + ) u+ - (27 - F»+) «+]^ 1/2 = S 3 (1) , 

3 ' (4.60) 

[u- (I- F< + ) 3* - (7 - F»+) «+] * +1 = S< 2 \ 

(4.61) 

[«+ (27 + F c+ ) u+ - (7 - 7™+) u+]”^ = 5 2 (2) , 

(4.62) 


[tT — (7 — F< + ) «+ + (27 + F v+ ) «+]'’* * = S 3 (2) , 

(4.63) 

where (;', fc,n -f 1/2) £ fli is assumed in Eqs. (4.58)- 
(4.60), while (i,t,n+l) £ ^ is assumed in Eqs. (4.61) 
-(4.63). Because 8 


(i) 


8 W 

s 2 » 


. 0 ) 

s 3 * 


s (2) .P) 


and « 3 2 \ de- 


note the expressions on the right sides of Eqs. (2.59)- 
(2.64), respectively, Eqs. (4.58)-(4.63) are the Euler 
images of the latter if one adds the following rule of 
substitution: 


Rule 4: 

be replaced by s£ k \ k = 1,2, and l = 
It follows from Eqs. (4.58)-(4.63) that (i) 

1,2,3. 


wr-ie 


(4.64) 

and 

(tT+) n+1/2 = - 
' 'i.* 3 


(4.65) 

where (j, k, n + 1/2) € fti; and (ii) 



«Wi 


(4.66) 

and 


s™ - §?') 

(4.67) 


where (j, jb,n-f 1) £ fi 2 - Eqs. (4.64)-(4.67) are the 
Euler images of Eqs. (2.67), (2.68), (2.70), and (2.71), 
respectively. Note that the Euler images of Eqs. (2.66) 
and (2.69) are also true [2]. These images are equiva- 
lent to Eqs. (4.54) and (4.55). Because the latter are 
easier to use, they will be used exclusively in the fol- 
lowing development. 

With the above preparations, an Euler solver cam 
now be defined. It consists of two marching steps. The 
first step is formed by Eqs. (4.54), (4.64), and (4.65); 
while the second is formed by Eqs. (4.55), (4.66), and 
(4.67). As explained earlier, k = 1,2, and 

£ = 1,2,3, become known after and tZ?* 1 are 

evaluated using Eqs. (4.54) and (4.55), respectively. 
This Euler solver has a two-way marching nature sim- 
ilar to that of the a scheme. As a result, it must be 
neutrally stable, (i.e., no numerical diffusion) if it is 
stable. Because it is reversible in time, this solver can- 
not model a physical problem that is irreversible in 
time, e.g., an inviscid flow problem involving shocks. 
Hereafter, this new Euler solver will be referred to as 
the 2D Euler a scheme. 


At this juncture, note that the 2D Euler a scheme 
is greatly simplified by the fact that tZ?* 1 / 2 and tZ?* 1 , 
respectively, can be directly evaluated in terms oi the 
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marching variables at the nth and (n + l/2)th time 
levels (see Eqs. (4.54) and (4.55)). As a result, the 

matrices and (E^*)?* 1 , which are non- 

linear functions of 1 ^ 2 and fi?* 1 , respectively, can 
be evaluated easily. In other words, nonlinearity of 
the above matrix functions does not cause a particu- 
lar problem for the Euler a scheme. 

To explain how Eqs. (4.54) and (4.55) arise, note 
that the conservation condtions (i) 

/ h' m • ds = 0 (4.68) 

J 5(C£?( 1 )(y,fc,n+l/2)) 

for any (j, Jb,n -f 1/2) £ fli, and (ii) 

f h* m -ds = 0 (4.69) 

J S(CJE?( a )(j ,Jb,n+l)) 


development, we have moved the center of SE(j, Jb, n-t- 
1/2) to the center of the top face, i.e., away from point 

G'. 

Next we shall construct the 2D Euler a-e scheme, 
i.e., the Euler version of the 2D a-e scheme. For this 
scheme, we shall use the CEs defined in Eqs. (3.1) and 
(3.2), i.e., Eqs. (4.68) and (4.69) will be assumed. Thus 
Eqs. (4.54) and (4.55) will also be part of the Euler a-e 
scheme. In the following, we shall describe the rest of 
the 2D Euler a-e scheme. 

As a preliminary, note that can be evalu- 

ated by a direct application of Eqs. (4.56) and (4.57), 
if one does not mind inverting the 4x4 matrices 

* and (^n^ + ) k • Alternatively, one 

may use the method of Gaussian elimination to ob- 
tain as the solution to 


for any (j, Jb,n 4* 1) € fi 2 , are the direct results of 
Eqs. (4.28) and (4.29), the basic assumptions of the 
Euler a scheme. According to Eq. (3.1), CE( l \j, Jb, n-f 
1/2) is the cylinder A! B' G Ft E' F' ABCDEF depicted 
in Fig. 7(a). Except for the top face A'B'GiyE'F 1 , 
the other boundaries of this cylinder are the sub- 
sets of three solution elements at the nth time level. 
Thus, for any m = 1,2, 3, 4, the flux of h* m leaving 
CE( l \j, Jb, n + 1/2) through all the boundaries except 
the top face can be evaluated in terms of the marching 
variables at the nth time level. On the other hand, 
because the top face is a subset of SE(,;, Jb,n -f 1/2), 
the flux leaving there is a function of the marching 
variables associated with SE(,;,Jb, n + 1/2). Further- 
more, because the outward normal to the top face has 
no spatial component, Eq. (4.25) implies that the to- 
tal flux of h leaving CE^^j, Jb, n -f 1/2) through the 
top face is the surface integration of over the top 
face. Because the center of SE(j , Jb, n -f 1/2) coincides 
with the center of the top face , it is easy to see that the 
first-order terms in Eqs. (4.15) do not contribute to the 
total flux leaving the top face. It follows that the to- 
tal flux leaving the top face is a function of (u m )?* 

only. As a result of the above considerations, 
can be determined in terms of the marching variables 
at the nth time level by using Eq. (4.68) only. Simi- 
larly, u?! 1 can be determined in terms of the march- 
ing variables at the (n + l/2)th time level by using 
Eq. (4.69) only. Eqs. (4.54) and (4.55) are the direct 
results of Eqs. (4.68) and (4.69), respectively. 

In an extension currently under development, the 
mesh used is not uniform in space. As a result, point 
G 1 depicted in Fig. 7(a) generally is not the center 
of the top face referred to earlier. To simplify the 


( 4 ; 


and 


( 2 ) 




(4.70) 


(4.71) 


t) 


To simplify computation further, in the following de- 


velopment of the Euler a-e 

scheme, we shall 

assume 

that (i) 





W'M* 

(o+r 
11 ) 

(4.72) 

where £ 

= 1,2,3 and (i,ib,n 

+ 1/2) 6 £2i, and 

(ii) 


on: = < 

' s (2)+\"+ 1/2 

(4.73) 

where £ 

= 1, 2, 3 and ( j , k, n 

+ 1) 6 ^2* 


To proceed, let (i) 



8 i — 


-(j+F> + )ff+] 

n 

O.M.i) 


[u + (2/-F c+ )tT ( H 


(4.74) 

r(l) 4£f 
S 2 ~ 

--(l+r> + )u + 

r 




(4.75) 

and 




s 3 — 

[tf- (/ + /’<+)« + 

+ ( 2 i - r> + ) u+ 

1" 

1 0\M,3) 




(4.76) 
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where (j, Jfc,n + 1/2) G fli; and (ii) 

«?’ =' [„-+ (/- jx+) «-+ + (/ - r> + ) “7]”* 1 ", 

'(4.77) 

sf> « [a - ( 2 / + F<+) « ( + + (I - F ’+) 8+] 

(4.78) 

and 

4 ” - [ ff+ <' - ^ + ) »'< + - < 2 ' + ^ ”‘- + ] „Z> 

(4.79) 

where (j, k, n + 1) 6 SV Then, with the aid of 
Eqs. (4.34)-(4.51), (4.72), and (4.73), Eqs. (4.56) and 
(4.57) imply that 

S^ k) = 4 k \ * = 1,2; £=1,2,3 (4.80) 

As a result, assuming Eqs. (4.72) and (4.73), Rule 4 
given above should be replaced by 

Rule 5: 8^ be replaced by s^ k \ k = 1,2, and l = 1,2,3. 
It follows that Eqs. (4.74)-(4.79) are the Euler images 
of the fact that s^, Jfc = 1,2, and £=1,2, 3, represent 
the expressions on the right sides of Eqs. (2.59)-(2.64). 

Combining Eqs. (4.64), (4.65) and (4.80), for any 
(i, k y n+ 1/2) G Oi, one has 



/ .\*+i/2 / 

' .\n+l/2 

(4.81a) 


("?),. =( 


and 

w>;r 

= (0,T' 2 

(4.815) 

where 


(.-<» - 4«) 

(4.82a) 

and 

/ 7 o+sn+l/2 dtf 1 
K ” 'it* ~ 3 

(4 ,, -4 ,) ) 

(4.826) 


Note that: (i) As a result of Eqs. (4.74)-{4.76), and 
(4.82a,b), ^ and (u®*)"* 1 ^ 2 can be evalu- 

ated in terms of the marching variables at the nth time 
level; and (ii) Eqs. (4.82a, b) are the Euler images of 
Eqs. (3.15a,b). 

Similarly, Eqs. (4.66), (4.67), and (4.80) imply 
that, for any (j, Jb,n + 1) £ ^ 2 , 



and 

W),? 


(4.836) 

where 

(s;+)" +1 * 1 

\ < )j,k 3 


(4.84a) 

and 

f.7°+'i n+1 d - I 1 
v " h,* ~ 3 1 

[.f’-«i 2) ) 

(4.846) 


Note that: (i) As a result of Eqs. (4.77)-(4.79) and 

(4.84a, b), k and 1 can evaluated 

in terms of the marching variables at the (n 4* l/2)th 
time level; and (ii) Eq. (4.84a,b) are the Euler images 
of Eqs. (3.27a, b). 

Furthermore, for any (j,k,n) G fi, let (t?t)"* 
denote the column matrix formed by ( u mt )"*, m = 
1,2, 3, 4. Then Eqs. (4.27), (4.30), and (4.31) imply 
that 


«« = - ^ (><+«+ + F* + 17+) . (4.85) 


With the above preliminaries, the Euler a-e 
scheme and the Euler weighted-average a-e scheme can 
be developed in a fashion parallel to the development 
of their non-Euler counterparts described in Sec. 3 [2]. 
Excluding a few exceptions to be discussed shortly, all 
other equations related to the Euler schemes are the 
trivial Euler images of its non-Euler counterparts un- 
der the substitution rules given in Rules 1-3, 5, and 
Rule 6: Any scalar variables be replaced by its column- 
matrix counterpart. 


In other words, one can obtain the Euler image from 
its non- Euler counterpart by adding an arrow over the 
symbols representing scalar variables. The complete 
set of equations that define the Euler a-e scheme and 
the Euler weighted-average a-e scheme is given in [2]. 
In the following, we shall only discuss those equations 
which cannot be obtained using the above substitution 
rules. One half of these special equations is formed by 
the Euler counterparts to Eqs. (3.62), (3.66a,b), and 
(3.67a, b). The other half is identical to the first half 
except that it is associated with (j, ib, n+ 1) G £^ 2 - As a 
result, we shall further restrict the discussion to those 
in the first half. 


Let the mth components of the 4x1 column 
matrices (fif*)# 1 '*, (u 1J W+ )^ 1/2 ) and 

(^)-t 1/2 , be denoted («$ + ® 1/2 . (^ + )£ 1/2 > 

(umr)?* 1 ^ 2 , and (timy)?* ^ 2 , respectively (these ma- 
trices can be defined by using the Euler images of 
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Eqs. (3.42M3.47), and (3.52)-(3.57)). Then, for 
m = 1,2, 3, 4, and £ = 1,2,3, the Euler versions of 
Eqs. (3.62), ( 3 . 66 a,b), and (3.67a, b) are (i) 

Omi = )/(«&) 2 + («$)* *, (4.86) 


(ii) 


= 0, if e ml -0m2 - e m3 = o, (4.87a) 

and 


u 


__ 

mC 


(e m 20 m3 ) a ^+(^ 3 ^ 1 )° +(g ml g m2 ) a tg- 

{OmlOml)* + (^m2^m 3 ) a + (Wml)° 

(4.876) 

otherwise; and (iii) 


its simplicity, particularly the fact that it does not use 
(i) any mesh refinement technique, or (ii) any mov- 
ing meshes, this new solver is capable of generating 
highly accurate unsteady shock solutions. The shock 
discontinuities are resolved almost to within one mesh 
interval. 

In this section, accuracy of the Euler weighted- 
average a-e scheme defined in Sec. 4 will be evaluated 
using a steady-state shock reflection problem [ 14]. The 
computation domain and the shock locations ( AE and 
EF) are depicted in Fig. 14. The lower boundary is a 
solid wadi. Assuming 7 = 1.4, the exact Euler solution 
to this problem is: 

(a) In the region ABE , 

u = 2.9, v = 0., p = 1.0, p = 1.0/1.4. 

(5.1) 

(b) In the region AEFD , 


* l rnr) — if ^ml — $m 2 — ^m3 — 0, (4.88a) 


and 


u 


v+ — 

mrj 


(^m2^m3) a (^m3^ml )** "K^ml ^2)° u ml) 

(Pml^m2) Q + (^m2^m3) a d* (^m3^ml) a 

(4.886) 

otherwise. Here (i) a > 0; (ii) the Euler counterpart 
to tij? + is the column matrix formed by m = 

1,2, 3, 4; and (iii) the Euler counterpart to ti^ + is the 
column matrix formed by m = 1 , 2 , 3, 4. 

This section is concluded with the following com- 
ments: 


(a) Because of the assumptions made in Eqs. (4.72) 
and (4.73), the Euler a scheme is not the special 
case of the Euler a-e scheme with e = 0. 

(b) Because (i) a fractional power is costly to evalu- 
ate, and (ii) evaluation of (9 m i) a does not involve 
a fractional power if a is an even integer, the Eu- 
ler weighted-average a-e scheme is more computa- 
tionally efficient if a is an even integer. 


5. Numerical Results 

In [11], several numerical solutions of Eq. (2.1) 
generated using the a-e schemes are compared with the 
exact solutions or the numerical solutions generated 
using traditional methods. These comparisons show 
that the a-e scheme, which includes the a scheme as 
a special case with e = 0 , is an accurate solver for 
Eq. ( 2 . 1 ). 

The a-e scheme was also generalized in [ 11 ] to 
solve the 2D inviscid Burgers’ equation. In spite of 


u = 2.6193, v = -0.50632, 

p = 1.7000, p = 1.5282. ( ’ 

(c) In the region ECF, 

u = 2.4015, v = 0., p = 2.6872, p = 2.9340. 

(5.3) 

Note that the Mach number is equal to (i) 2.9 in the 
region ABE ; (ii) 2.3781 in the region AEFD] and (iii) 
1.9424 in the region ECF. 

The mesh used in the current numerical calcula- 
tions is depicted in Fig. 15. Again a mesh point £ fii 
is marked by a solid circle; while a mesh point £ &2 
is marked by am open circle. The mesh is a speciad 
case of that depicted in Figs. 3-6 with 6=0. Note 
that (i) only the mesh points £ Q 2 are present at the 
inflow boundary, and (ii) the mesh parameter w is so 
chosen that only the mesh points £ Q 2 are present 
at the outflow boundary. Moreover, for simplicity, a 
mesh point and the corresponding maLrching variable 
will be identified by the time-level number n, amd two 
new mesh indices r and s which are given in Fig. 15 as 
a pair of integers enclosed in parentheses. Note that, 
for the mesh points £ fix, r = 1,2,3, ... ,iZ, iZ + 1, 
and 8 = 1,2,3,..., 5. On the other hand, for the 
mesh points £ fi 2 , r = 1,2,3, + 1, and 

sal, 2, 3, ...,5,5+1. Obviously two different mesh 
points at the same time level adways have different 
padre of r and s. 

In the current numerical calculation, at the time 
level n = 0, ti m , m = 1,2, 3, 4, at all mesh points aire 
calculated using Eq. (5.1). Also we assume that 

= 0. m — 1,2, 3,4. (5.4) 
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The above initial conditions are also assumed at the 

inflow boundary for all n = 1,2,3, At the upper 

boundary, for all n = 1/2, 1,3/2, 2, . . Eq. (5.4) is also 
assumed. Moreover, Um, m = 1,2, 3, 4, are calculated 
using Eq. (5.2). 

To impose the proper boundary conditions at the 
lower bou ndar y, note that the solid wall boundary con- 
ditions at BC (see Fig. 14) are equivalent to the condi- 
tion that the flow field below BC is the mirror image 
of that above. By using Eq. (4.1) and the fact that 
y = 0 at any point on BC y it can be shown that the 
last condition implies that 


u 3 (x,-y) = -u 3 (x,y) 
du m (x ,-y) _ drimix^) 


dx dx ’ 

du 3 (x,-y) du 3 (x,y) 


dx 


dx 


dy 


dy 


and 


du 3 (x,-y) du 3 (x,y) 


dy dy 

Consider the mesh depicted in Fig. 15. Then it be- 
comes clear that the numerical analogues of Eqs. (5.5a) 
-(5.7b) are 

( Um )fl+i,» = m = (5.8a) 


(“Wkm,. = (“£<)*, • m = I - 2 ' 4 (5- 12 «) 

and 

R)k + 1i , = - (4) r ., 

respectively. Equations (5.8a, b), (5.11a,b), and 
(5.12a,b) are the boundary conditions at the lower wall 
(a solid wadi) in the current numerical calculations. In 
other words, the marching variables associated with 
the mesh points below the solid wall will be deter- 


(5.5a) 

mined using these equations. 


Next we discuss the outflow boundary conditions. 

(5.56) 

For any n = 1,2,3,..., and r = 1,2,3, ..., R, 
sume that 

we as- 

(5.6a) 

(“ m )r,s+l = ( Um )r,S > m= 1> 2 >3> 4 

(5.13) 

(5.66) 

(«mx)" 5+ i = °, m = l,2,3,4 

(5.14) 

(5.7a) 

and 



(«my) r 5 +1 = ( u my) rS » m = 1,2, 3, 4 

(5.15) 

(5.76) 

When the time-marching solution reaches its steady- 


state limit, the above conditions can be considered as 
a result of the requirement that the partial derivatives 
of the flow variables with respect to x are zero at the 
outflow boundary. Because b = 0, it can be shown 
that Eqs. (5.14) and (5.15) are equivalent to 


(5.86) 

(“mx)* + 1,. = M)l,.< m = 1 > 2 > 4 ( 5 - 9a ) 

Ml+i,. = ~Ml,. ( 5 - 9fe ) 

( u my) fl+l t = ~( u my) R t , m = 1, 2,4 (5.10a) 

and 

Ml +1 , = Ml,. (5-106) 

respectively. According to Fig. 15, the range of s in 
Eqs. (5.8a)-(5.10b) is dependent on the time level n. 
Let (i) 5 + = f S -f 1, and S~ = f S if S is even; and 

(ii) S+ d = S, and S~ = f S - 1 if S is odd. Then (i) 
8 = 2,4, 6,. ..,5“ if n = 1/2, 3/2,..., and (ii) 8 = 

1,3,5,..., S + if n = 1,2, Furthermore, because 

6 = 0, it can be shown that Eqs. (5.9a,b) and (5.10a, b) 
are eqivalent to 

(«mc)* +1> , = (Ml,. > rn = 1,2,4 (5.11 a) 


and 




Wo)" 


5+1 


= \ (“»( 

= j 



n— 1/2 
r,5 



n— 1/2 
r »5 


(5.16) 


(5.17) 


where m — 1,2, 3,4, n = 1,2,3,..., and r = 
l,2,3,...,iJ. Equations (5.13), (5.16), and (5.17) are 
the outflow boundary conditions in the current numer- 
ical calculations. As a result, the marching variables at 
the outflow boundary will be determined using these 
equations. 

With the aid of the above initial and boundary 
conditions, the marching variables at all time levels 
can be determined using the Euler weighted-averaged 
a-e scheme. As an example, at any n = 1/2, 3/2, . . ., 
the marching variables associated with the mesh point 
(2, 1) (marked by a solid circle in Fig. 15) can be de- 
termined in terms of those associated with the mesh 
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points (1, 1), (2, 1), and (2,2) at the (n — l/2)th time 
level (marked by open circles). As smother example, 
at any n = 1,2,3,..., the marching variables associ- 
ated with the mesh point (1,3) (marked by an open 
circle) can be determined in terms of those associated 
with the mesh points (1,2), (2,3), and (1,3) at the 
(n — l/2)th time level (marked by solid circles). 

According to Fig. 14, the distance between the 
inflow and the outflow boundaries is 4., while the dis- 
tance between the upper and the lower boundaries is 
1.. On the other hand, according to Fig. 15, the above 
two distances are w • S and 2 h • fZ, respectively. Thus 

W =T}> and h= jR ( 518 ) 

In addition to the initial conditions, the boundary 
conditions, and the integers R and 5, the other input 
parameters for the current numerical calculations are 
e, a, Af, and a positive integer n t . Here we assume 
that the time marching ends at the n t th time level, 
i.e., at t = T = f n t • At. 

It is shown in [2] that, for any Euler solver con- 
structed in Sec. 4, a local CFL number v t associated 
with any mesh point (j, Jb, n) £ fl can be defined in 
terms of tx, v, c, tn, h, and At. Here tx, v, and c are the 
r-velocity, the y- velocity, and the sonic speed at the 
mesh point, respectively. Two global CFL numbers 
are considered in the current calculations. The first, 
denoted by v em „ is the maximum of u t with respect 
to the steady-state solution given in Eqs. (5.1)-(5.3). 
The second, denoted by i/ emi is the largest value of v t 
ever reached at any mesh point (j, Jfc,n) £ fi, where 
n = 0, 1/2, 2, 3/2,2,..., n t . Excluding the initial and 
the boundary conditions, u em9 is dependent on R, 5, 
and At only. On the other hand, v em is a function of 
R y Sy At, n t y Cy and a. According to a series of nu- 
merical experiments, the value of a plays only a minor 
role on the stability of the Euler weighted-average a-e 
scheme. Generally the scheme is stable if 

0 < e < 1 and i/ em < 1. (5.19) 

To measure the convergence of a time-marching 
solution to the corresponding steady-state solution 
(note: this steady-state solution generally differs from 
the exact solution given in Eqs. (5.1)-(5.3)), for any 
n = 1,2, 3, ...,n t , and m = 1,2, 3,4, let 


Here, for any m = 1,2, 3, 4, c™ is the maximal value 
of It^l within the exact steady-state solution defined 
by Eqs. (5.1)-(5.3). It can be shown that c\ = 2.6872, 
ci = 6.4534, C 3 = 0.86073, and C 4 = 15.084. Moreover, 



if s is odd; 
otherwise. 


(5.21) 


According to Fig. 15, the summation that takes place 
in Eq. (5.20) involves all the mesh points at the nth 
time level excluding those located (i) at the inflow 
boundary, (ii) at the upper boundary, and (iii) below 
the lower boundary. Because n = 1,2,3, ... ,n*, the 
mesh points involved in the summation are all marked 
by open circles in Fig. 15. The values of tx m at the 
inflow and the upper boundaries do not change with 
time, while those at the mesh points in (iii) are depen- 
dent on the values of tx m at other interior mesh points. 
Note that the values of tx m at the outflow boundary 
change with time and are dependent on those at a 
lower time level. Because the summation involves a 
total of R x S mesh points, the result of this summa- 
tion divided by Rx S is the average value of the change 
of tx m at the same spatial mesh point (measured by the 
absolute value of this change) from the (n — l)th time 
level to the nth time level per mesh point. This av- 
erage value is further normalized using the constant 
Cm. If we further assume that the time-marching solu- 
tion converges to a steady-state solution that is similar 
to the exact steady-state solution (such that the nor- 
malization by c m makes sense), then E m (n) can be 
interpreted as the average number of correct signifi- 
cant figures in tx m at the nth time level as compared 
with the converged value of tx m (which, of course, is 
not identical to that given in Eqs. (5.1)-(5.3)). 

Because the time marching solution can not reach 
a steady-state solution before the boundary conditions 
are fully felt at sill interior points, rapid convergence 
generally can not occur before the time has elapsed 
that allows a fluid particle to travel the full length 
of the computation domain. It can be shown that, 
for the solution given in Eqs. (5.1)-(5.3), the average 
value of tx over the computational domain is 2.6261. 
Thus an average fluid particle requires 4.0/2.6261 = 
1.5232 time units to travel from the inflow boundary 
to the outflow boundary. The number of time steps 
corresponding to the above number of time uints is 


E m {n )® 


( 1 

’ 5+1 

r(,)+R-l T | 

E 



1=2 

r = r (*) JJ 


(5.20) 


def 1.5232 

n c = 

c a t 

i.e., rapid convergence can not occur before n > n c . 

With the above preliminaries, the numerical re- 
sults generated using the Euler weighted-average a-e 
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scheme can now be presented. Six test problems, with 
different combinations of c, a, R , 5, a i, and n*, are 
defined in Table 1. For each problem, the values of T, 
i/ em „ i'em, and n c are also given in the same table. In 
Figs. 16-21, the numerical results (triangular symbols) 
of the pressure coefficient c p at n = n t for Problems 
#l-#6 are compared with the exact solution (solid 
lines). Here 

(5 - 23) 

with Moo = 2.9 and Poo — 1.0/1.4 being the in- 
flow Mach number and pressure, respectively. Note 
that: (i) at the mid-section of the computation do- 
main (y = 0.5 in Fig. 14), two neighboring mesh points 
at the same time level are separated by a distance 
= 2ty, and (ii) the mesh points at the n t th time level 
are marked by open circles in Fig. 15 because n t is a 
whole number. In Figs. 16-21, the values of E m (n), 
m — 1,2, 3, 4, are also plotted against n for all six 
test problems. In Fig. 22, twenty-six pressure con- 
tour levels between the values of 0.6 and 3.1 with uni- 
form increment 0.1 were used for the contour plots of 
Problem #3. Finally, for Problem #3, a 3D pressure- 
distribution plot is shown in Fig. 1. Note that, for 
easier visualization, only pressure values at the mesh 
points with 8 = 2,4, 6, . . . are used in the last plot. 

The significance of the results shown in Figs. 1 
and 16-22 is discussed in the following remarks: 

(a) From Table 1 and the results shown in Figs. 16-18, 
it appears that the convergence to steady-state is 
much faster with a smaller value of v tm (or i/ cm# ). 
As a matter of fact, convergence to steady-state 
can reach a plateau representing some number of 
correct significant figures if v tm is too close to 1. 
From Table 1 and a comparison among Figs. 16, 
19, and 20, one also concludes that slower con- 
vergence generally occurs with a value of c much 
smaller than 0.5. A comparison between Figs. 16 
and 21 reveals that a change of the value of a 
from 2 to 1 also causes a slight decrease in con- 
vergence rate. Because numerical diffusion gen- 
erally increases with (i) a smaller value of i/ cm , 
(ii) a larger value of e, and (iii) a larger value 
of a, one may conclude that faster convergence 
generally occurs with larger numerical diffusion. 
This trend is consistent with the fact that shocks 
cannot be formed without physical or numerical 
diffusion. 

(b) The effectiveness of weighted- averaging as a tool 
to surpress numerical oscillations near discontinu- 
ities is clearly demonstrated by the results shown 


in Figs. 1 and 16-22. Moreover, the current 
weighted-averaging does not cause the smearing 
of shock discontinuities and has no discernible ef- 
fect on the smooth part of the solution. From 
table 1 and a comparison between Fig. 16 and 21, 
one also concludes that the increase of the value 
of a from 1 to 2 has a marginal impact on the 
numerical results. 

(c) Comparing the numerical results shown in 
Figs. 16-21 with the exact solution, one concludes 
that the Euler weighted-average a-c scheme is ca- 
pable of generating highly accurate solutions for 
the steady-state shock reflection problem under 
consideration. Also a comparison of the results 
shown in Figs. 16, 17, and 19-21 reveals that ac- 
curacy of the numerical results generally is not 
sensitive to the change of the values of v tm , e, 
and a. An exception is that numerical results 
may become more diffusive and thus shock reso- 
lution becomes less sharp if the value of e is too 
large, e.g., e == 0.8 in Problem #5. Finally, a com- 
parison of the results shown in Fig. 18 (Problem 
#3) with the results of other test problems reveals 
that accuracy increases sharply with a decrease of 
the mesh size. 

This section is concluded with a brief discussion 
on recent applications of the current solver to com- 
putational aeroacoustics (CAA) problems. CAA ik an 
area of current interest in CFD. It is of both theoretical 
interests and practical importance. For a scheme to be 
a useful research tool, it must be accurate enough to 
resolve sound wave details. Furthermore, its boundary 
conditions must be non-reflecting. There have been a 
large number of papers published on these two topics. 
In practice, CAA will help to reduce jet noise level 
for air-borne vehicles, and hence becomes an impor- 
tant topic in HSR (High Speed Research) and AST 
(Advanced Subsonic Technology) programs. 

The most popular numerical schemes for CAA are 
the high order (4-6th) compact or non-compact dif- 
ference schemes, marching forward by Runge-Kutta 
method. These schemes work quite well when incorpo- 
rated with delicatedly designed non-reflecting bound- 
ary conditions. They are even capable of capturing 
Mach waves and weak shocks. 

When conducting numerical experiments of the 
current solver for CAA problems, surprisingly we 
found the following attractive features: 

• Accuracy of the current solver is comparable to 
that of a 6th-order compact difference scheme, 
even though nominally the current solver is only 
of 2nd-order accuracy. For example, the smallest 
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eddy can be resolved in 3-4 grid cells. 

• Generally, the non-reflecting (radiation) bound- 
ary condition can be implemented in a simple way 
without involving characteristic variables. 

• Most importantly, the current solver is capable 
of handling BOTH continuous and discontinuous 
flows very well and thus provides a unique numer- 
ical tool for solving those flow problems where the 
interactions between sound waves and shocks are 
important, such as the noise field around a super- 
sonic over- or under expansion jet. 

Details of our investigation will be reported else- 
where [13]. Here as an example, we present the nu- 
merical results of a calculation involving a free shear 
layer which is subjected to upstream time-dependent 
perturbation. The Mach number of the fast stream is 
1.5 and the slow stream is subsonic. Fig. 23 illustrates 
the contours of pressure, v-velocity and vorticity in 
the near field. The computational mesh is formed by 
240 x 150 uniform cells. From v-velocity contours, ed- 
dies created by the upstream perturbation are clearly 
visible. On the subsonic side, sound waves propagate 
freely to far field, while on the fast stream side (lower 
1/3 domain) the sound waves are bounded by the Mach 
line, even smaller eddies are created around the Mach 
line due to their interactions. In summary, the current 
results reveals a lot of physical details of the perturbed 
shear layer and we believe that the current method will 
develop into a robust and unique numerical technique 
for aeroacoustics computation. 

6. Conclusions and Discussions 

A new numerical method is being developed for 
solving one-dimensional and multidimensional flow 
problems. This new method represents a clear break 
from the traditional methods in the basic concept 
of numerical discretization . It emphasizes simplicity , 
generality and accuracy . The history of this new 
method and the considerations that motivate its de- 
velopment are described in Sec. 1. 

In this paper, the same design principles that were 
used to construct several solvers for ID time-marching 
problems [5] are used to construct their 2D counter- 
parts. Because of the similarity in their designs, each 
of the present 2D solvers shares with its ID counter- 
part virtually the same fundamental characteristics. 
Furthermore, it has been shown that the 2D solvers, 
as in the case of the ID solvers, generally are more ac- 
curate than the traditional solvers despite the advan- 
tage the current solvers have over the latter in simplic- 
ity and generality. Accuracy of the current 2D Euler 


solver is most vividly demonstrated by the pressure- 
contour plot (Fig. 22) and the 3D pressure-distribution 
plot (Fig. 1) it generates for a famous shock reflec- 
tion problem [14]. Both the incident and the reflected 
shocks are resolved by a single data point without the 
presence of numerical oscillations near the discontinu- 
ity. 

Construction of the ID solvers referred to above 
is simplified by the use of a mesh that is staggered in 
time [1,9]. Its use results in the simplest stencil possi- 
ble, i.e., a triangle in 2D space-time with one vertex at 
the upper time level and other two at the lower time 
level. Similarly, construction of the current 2D solvers 
is simplified by the use of a nontraditional space-time 
mesh that is also staggered in time (Figs. 3-6). Its use 
results in the simplest stencil possible, i.e., a tetrahe- 
dron (Fig. 10) in 3D space-time with one vertex at the 
upper time level and the other three at the lower time 
levels. 

The meshes used by the ID and 2D solvers con- 
sist of whole-integer and half-integer time levels with a 
half-integer time level being sandwiched between two 
whole-integer time levels, and vice versa. The spa- 
tial positions of the mesh points at a whole-integer 
(half-integer) time level coincide with those at an- 
other whole-integer (half-integer) time level. However, 
the spatial positions of the mesh points at a whole- 
integer time level shift from those at a half-integer 
time level. For the mesh used by the ID solvers, the 
spatial projection of a mesh point at a whole-integer 
time level is right at the center of those of two neigh- 
boring mesh points at a half-integer time level, and 
vice versa [1,5,9]. It follows that the stencil of the ID 
solvers is always an isosceles triangle, i.e., one cannot 
distinguish a stencil with its upper vertex at a whole- 
integer time level from another with its upper vertex 
at a half-integer time level. As a result, each of the ID 
solvere constructed in [1,5,9] is formed by two iden- 
tical marching steps. Contrarily, for the present 2D 
solvere, a stencil (a tetrahedron) with its vertex at a 
whole-integer time level is different from another with 
its vertex at a half-integer time level (Fig. 10). Thus, 
each of the present 2D solvers is formed by two dis- 
tinctly different marching steps. In spite of their struc- 
tural differences, the last two marching steps compen- 
sate each other and their combination results in several 
important symmetric properties that are discussed in 
[ 2 ]- 

The Euler a scheme constructed in Sec. 4 is free 
from numerical diffusion when it is stable. This scheme 
is a limiting case of a Navier-Stokes solver currently 
under development, i.e., the former is a special case of 
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the latter when the viscosity vanishes. As a result, the 
new Navier-Stokes solver will have a special property 
that a classical solver lacks, i.e., as the physical diffu- 
sion (viscosity) approaches zero, so does the numerical 
diffusion. The significance of this property was dis- 
cussed earlier. Because a Navier-Stokes problem fun- 
damentally is an initial-value/boundary-value prob- 
lem, i.e., information from any spatial point can be 
felt instantly by other spatial points, the new Navier- 
Stokes solver is implicit when viscosity is present. 
However, it becomes explicit when viscosity is absent. 

Finally, note that a new implicit solver for 
Eq. (1.3) has been developed recently [4]. This new 
solver shares with the above Navier-Stokes solver es- 
sentially the same characteristics. In the inviscid limit, 
this new scheme becomes explicit and its amplification 
factors are identical to those of the Leapfrog schemes. 
On the other hand , in the pure diffusion limit , its prin- 
cipal amplification factor becomes the amplification 
factor of the Crank-Nicolson scheme. By using an ap- 
proach similar to that described in [4], a new a-e solver 
for Eq. (1.1) was also developed. Stability of this new 
a-c scheme is again limited by the CFL condition and 
0 < e < 1. Moreover, if c = 0, its amplification factors 
are identical to those of the Leapfrog scheme. On the 
other hand , if c = 1, its principal amplification factor 
becomes the amplification factor of the Lax-Wendroff 
scheme . This new a-c scheme will be reported in the 
near future. 
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Figure 2 .— The leapfrog scheme using a moving mesh and 
the upwind leapfrog scheme. 



Figure 3. — The relative spatial positions of the mesh points 
and the mesh points €D 2 (dash lines are spatial 
boundaries of the conservation elements depicted in 
figs 7(a) and 7(b)). 



Figure 4. — The spatial mesh indices (j, k) of the mesh points 
6^ (n = ±1/2, ±3/2, ±5/2, •••). 



Figure 5. — The spatial mesh indices (j, k) of the mesh points 

€O,(n=0,±1,±2,-). 
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Figure 6. — The spatial mesh positions of the mesh points 
marked by • and those marked byo . 



G' = G» k, n + i) 

k+ |- n+ f) 

D=(j -i, k .i >n+ l ) 

F = 0 + | k -l.n + I) 

CE^ 0, k, n + 1) - box CDEGC'D'E'G' 
CE^ 0. k. n 1) - box AGEFA'G'ET' 
CE^ fl, k, n + 1) - box ABCGA'B'C'G' 


G' = G. k, n + 1) 

A' = G+ 1, k + l,n + l) 
3 3 

|,k+ 1, n + 1) 
E '=fl+ |.k- n + 1) 


SE (2) (j, k, n ♦ 1) = the union of four planes 
A'B'C'D'ET', GG "A "A, GCC' G", and GG"E"E 
and their immediate neighborhoods. 

Figure 8. — (a) Conservation elements CE ( p (j, k, n + 1), € = 1, 
2, 3, j, k, = 1/3, 1/3 ±1,1/3 ±2, and n = 0, ±1 , ±2, •••. 

(b) Solution elements SE^ (j, k, n + 1), j, k = 1/3, 1/3 ±1, 

1/3 ± 2 . — , and n = 0, ±1 , ±2, •••). 



Figure 7. — (a) Conservation elements CE ( P (j, k, n + 1/2), 

€ = 1 , 2, 3, and j, k, n = 0, ±1 , ±2, — . (b) Solution elements 
SE (1) 0, k, n + 1/2), j, k, n = 0, ±1, +2, •••). 





Figure 9. — Geometry of the hexagon ABCDEF. (a) Relative 
positions of the vertices in terms of (x, y). (b) Relative 
positions of the vertices in terms of (j, k). (c) Relative 
positions of the vertices in terms of (£, x|). 
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n), (j - 2/3, k + 1/3, n) and (j + 1/3, k - 2/3, n) with 0, k, n + 
1/2) (b) The mesh points (j, k, n + 1), (j - 1/3, k - 1/3, 

n +1/2), 0 + 2/3, k - 1/3, n + 1/2), and 0 - 1/3, k + 2/3, n + 
1/2) with 0, k,n + 1)eft 2 . 



Figure 12. — The space, (a) 0, k, n + 1/2) eH-j. (b) (j, k, 

n + 1) €H>2 a 




' n+ l 
J 3’ K + 3 


n + 1 


Figure 1 3. — The weighted-average a-€ scheme, (a) (j, k, n + 
1/2) € H-j. (b) Q, k, n + 1) € fl 2 . 
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Figure 1 4. — The computation domain and the shock 
locations of a steady-state shock reflection problem. 
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Figure 1 5. — ^The spatial locations and the new mesh indices 
(r, s) of mesh points (R = S = 4). 
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Figure 17. — Numerical results and convergence histories 
for problem #2. (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). (b) Converg- 
ence histories for u m , m = 1 , 2, 3, 4. 




Figure 16. — Numerical results and convergence histories 
for problem #1 . (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). (b) Conver- 
gence histories for u m , m = 1 , 2, 3, 4. 




n 


Figure 1 8. — Numerical results and convergence histories 
for problem #3. (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). (b) Converg- 
ence histories for u m , m = 1 , 2, 3, 4. 
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Figure 1 9. — Numerical results and convergence histories 
for problem #4. (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). 

(b) Convergence histories for u m , m = 1 , 2, 3, 4. 


Figure 21. — Numerical results and convergence histories 
for problem #6. (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). 

(b) Convergence histories for u m , m = 1, 2, 3, 4. 




Figure 20. — Numerical results and convergence histories 
for problem #5. (a) Pressure coefficients at the mid-section 
of the computation domain (y = 0.5 in Fig. 14). 

(b) Convergence histories for u m , m = 1 , 2, 3, 4. 



X 

Figure 22. — Pressure contours for problem #3. 
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Figure 23. — Aeroacoustic computation of a free shear 
layer. (The fast stream lies in the bottom 1/3 domain 
while the slower stream lies in the top 2/3 domain. 

A small sinusoidal time-dependent perturbation is 
applied at the shear layer separating two streams.) 
(a) Pressure contours, (b) v-velocity contours. 

(c) Vorticity contours. 


Table 1. — Definitions of test problems numbers 1 to 6 and the 
corresponding values of T, v^, and n c 



e 

a 

R 

S 

'At 

A 

* T 

V em 

V «n 

°c 

1 

0.5 

2 

20 

60 

0.01 

600 

6 

0.585 

0.6204 

152.32 

2 

.5 

2 

20 

60 

.015 

400 

6 

.8775 

.9305 

101.55 

3 

S 

2 

40 

120 

.0075 

800 

6 

.8775 

.9302 

203.09 

4 

2 

2 

20 

60 

.01 

600 

6' 

.585 

.6303 

152.32 

5 

.8 

2 

20 

60 

.01 

600 

6 

.585 

.6212 

15232 

6 

5 

1 

20 

60 

.01 

600 

6 

.585 

.6206 

15232 
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